{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d6d92d2c",
   "metadata": {},
   "source": [
    "# 习题\n",
    "## 习题19.1\n",
    "![image.png](./images/exercise1.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "09ac9a61",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "抽样样本数量: 1000000.0\n",
      "近似解: 2.5092124834947285\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "class MonteCarloIntegration:\n",
    "    def __init__(self, func_f, func_p):\n",
    "        # 所求期望的函数\n",
    "        self.func_f = func_f\n",
    "        # 抽样分布的概率密度函数\n",
    "        self.func_p = func_p\n",
    "\n",
    "    def solve(self, num_samples):\n",
    "        \"\"\"\n",
    "        蒙特卡罗积分法\n",
    "        :param num_samples: 抽样样本数量\n",
    "        :return: 样本的函数均值\n",
    "        \"\"\"\n",
    "        samples = self.func_p(num_samples)\n",
    "        vfunc_f = lambda x: self.func_f(x)\n",
    "        # np.vectorize 将一个对单个元素操作的函数转换为对整个数组或多个数组进行元素级操作的函数。\n",
    "        vfunc_f = np.vectorize(vfunc_f)\n",
    "        y = vfunc_f(samples)\n",
    "        return np.sum(y) / num_samples\n",
    "\n",
    "def func_f(x):\n",
    "    \"\"\"定义函数f\"\"\"\n",
    "    return x ** 2 * np.sqrt(2 * np.pi)\n",
    "\n",
    "\n",
    "def func_p(n):\n",
    "    \"\"\"定义在分布上随机抽样的函数g\"\"\"\n",
    "    return np.random.standard_normal(int(n))\n",
    "\n",
    "\n",
    "# 设置样本数量\n",
    "num_samples = 1e6\n",
    "\n",
    "# 使用蒙特卡罗积分法进行求解\n",
    "monte_carlo_integration = MonteCarloIntegration(func_f, func_p)\n",
    "result = monte_carlo_integration.solve(num_samples)\n",
    "print(\"抽样样本数量:\", num_samples)\n",
    "print(\"近似解:\", result)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0d42c9b0",
   "metadata": {},
   "source": [
    "## 习题19.7\n",
    "![image.png](./images/exercise7.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a43af0c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "均值: 0.4126647495883334\n",
      "方差: 0.018624023882369674\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi0AAAGvCAYAAACXeeU8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuO0lEQVR4nO3de1RU9f7/8dfgEAwIAyoGGIqKlZ2vBpnYKThqeSnTMu2ipeZKCy2rr3a1UhDt6LHSah39dsqUPKUl3TRPxy6rtDJLFNQu5jHFCyqXVKBQue7fH/6cE6nA4HD5wPOx1qzFns/e+/Oe+YD75d6f2WOzLMsSAABAI+fV0AUAAADUBKEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQXwgJSUFNlsNt1www2u55599lnZbDaNHTu24QqrQ3v27JHNZtOePXvOuk5kZKRSUlI83ndSUpJsNptOnDjhem7t2rWy2WzasmWLx/urztixY5vsOAONib2hCwCakq1bt7p+3rZtW6328f7770uShg4d6oGKGtZ7772nCy64oKHLOCc1GY9p06bVTzFAM0doATxo3759ys/PV1BQUKUA446mFFpiYmIauoRzVpPx6Ny5c/0UAzRzXB4CPKRDhw4KCwvTtm3bVFpaqp9++km9evVq6LIAoMkgtAAe1K1bN23btk3bt2+XzWZTly5dKrW/88476tatmxwOhy677DJ9/vnnrrbIyEjZbDa99tpreu2112Sz2WSz2bR27VrXOikpKQoKClJpaammT5+ujh076vnnn6/Ux+uvv64LL7xQPj4+io2N1ZdffulqS0pKUq9evTRlyhQFBgaqQ4cOmjNnjn7/Ze9lZWV64oknFBoaKn9/f910003av39/rd6Ps81pKS4u1v3336/w8HD5+/srPj6+VnNRysvLVVZWprKyMpWXl1dqKyoq0oQJE3T++efL6XTq2muv1e7duyut89JLL6lLly5yOBzq2rWrVqxYUan26sbjlLPNaenTp4+SkpK0ePFiRUZGKjAwULfffnuluTj79u3TwIEDFRAQoF69emnevHm68MILNWvWLI++V0BTQGgBPKhbt27aunWrtm3bpksuuUQtWrRwta1du1a33HKLhgwZon/961/q0aOHrr32Wu3YsUPSyfkfGzZs0KBBgzRo0CBt2LBBGzZs0GWXXXZaP3feeafWrFmjBx54QH369HE9v3TpUo0ZM0ZDhw7VqlWrFBkZqX79+mnz5s2udTZt2qRNmzbprbfe0oQJE/TUU09p4cKFrvaEhAQtWLBAiYmJWr58uXbv3q3evXursLDQY+/T7NmztXjxYs2ZM0dvv/22WrZsqdtuu83t/bRs2VLe3t7y9vZWv379KrVNnjxZb775pv7+97/r7bffVlFRke6++25X+9q1azVx4kSNHDlSq1ev1nXXXac77rhD+/btk+TeeFRl5cqV+utf/6pnnnlGc+fO1YoVK/TKK6+42seNGyebzaaVK1fqoosu0syZM5WSkqJbb73Vo+8V0BQwpwXwoG7dumnhwoUKCgpS9+7dK7XNmDFD119/vZKTkyVJcXFxev/997V8+XIlJSW55n+EhIRIkq644ooz9lFQUKBffvlF69evl7e3d6W26dOn64477tDcuXMlSf3799ell16qWbNm6b333pMk2e12rVixQqGhobruuuv0448/6oUXXtB9992nzMxMLVmyRK+88orGjRsn6eS8lC5dumjJkiV68MEHPfI+7dmzR5GRkRozZowkKTo6Wunp6bIsSzabrcb7Wb9+vc477zxJ0ubNmzVhwgRXW9++fTVy5Ej17dtXkrR9+3Y9/vjjlWqQpEmTJqlt27aKj4/X1VdfLT8/P0mq8XhUZ8eOHfrpp5/Uvn17SdIHH3xQaZL2hg0blJqaqquvvlphYWH65z//qY4dOyosLMxVpyfeK6ApILQAHtStWzd9//33CgwM1LXXXqvvvvvO1bZt2zYdOXLktKCxc+dOt/rw8vLSCy+8cNp+8vLytHfvXiUlJVVat2/fvnr33Xddz3Xu3FmhoaGu5djYWL355psqKyvTpk2bZFmWrrnmGld7RESELrzwQqWlpblVZ1XGjx+vt99+WzExMerTp4/i4+M1ZMgQtw/Cl112mXx9fSVJv/32W6W2m2++WcuWLdOYMWO0YcMG7dq1q9JlsBtuuEEdO3ZUjx49NHDgQF1xxRUaNmyYWrVqde4v8HeGDh3qCizSyRBUWlrqWr744ou1evVqXXXVVVq5cqWCg4N1/vnnu9o99V4BTQGXhwAPuuSSS1RcXKxPP/1U3bp1O619woQJSktLq/Q4dealpgICAtS1a9fTnv/9Afn3bDZbpbY/rldRUeGar1HTfZyruLg47dq1Sw899JCOHz+ucePGqXfv3iorK/PI/svLy3X11VcrKSlJnTt31vPPP6+VK1dWWqdVq1b64YcftHDhQrVu3Vpz587VRRdd5Lo85CnVfbIoJiZGKSkpcjqdmj17tlJSUuTl9d9/muv6vQJMQmgBPMjX11dRUVGSdNrlof/5n/9Rbm6uLr/8ctdj9erVWrNmzWn7OH78uNt9t23bVhEREfrss89cz1VUVOjzzz/X5Zdf7nru559/1oEDB1zLGzZsUOfOndWiRQv16NFDkirtIysrSzt27Ki0j3M1e/Zs/fTTTxo1apReeuklLV++XBs2bKh0ZupcfP/99/rqq6/06quvKjExUddff72ysrIqrfP222/rjTfe0JAhQ/S3v/1N6enpys/P1zvvvFNpvdqOxym/n9f0RxkZGVq1apUOHz6sH3/8UYcOHap0g0Kp7t8rwCRcHgI8rFu3bjp8+HClSzDSyfkm/fv319SpU9W/f3998803Sk5OVmpqaqX1YmNj9dhjj+lf//qXJGn//v2V5mpUZcaMGRo/frzatWunvn37avHixfrpp5+0aNGiSuvdcsstmjZtmtLT05WamuqaiNu5c2eNGTNGDz30kEpKShQeHq7ExESFhobqrrvuqu1bcprMzEy98cYbSk5OVqtWrbRkyRJ5e3urXbt2Htl/cHCwbDabli1bJpvNpk8//VTPPfecpJOfjrLb7SopKdGUKVMkSRdeeKG++OILlZWVqVOnTpX2dS7jUR2Hw6HDhw9r4cKF6tGjh44dO6awsDCFh4e71qnr9wowigXgnC1ZssTq0KGDZVmWNWPGDKtv376WZVnWnXfead15552u9d566y3rT3/6k+Xj42NdfPHF1pIlS07bV1lZmTVp0iQrODjYcjgcVkJCQqV+nE5nlbWkpKRYUVFRlre3t3X55Zdba9eudbUlJiZavXr1sqZOnWoFBgZabdq0saZPn25VVFS41ikpKbEeffRRKyQkxHI4HNYNN9xg7dmz57R+MjMzLUlWZmbmWWvp0KHDGV9jQUGBdc8991hhYWGWj4+P1b17d2vlypVVvq7fS0xMtCRZx48fdz33+eefW5KsjIwMy7Is6+WXX7bat29v+fj4WH/5y1+s1157zZJkffbZZ65t5s6da3Xp0sXy9fW12rdvbz399NOn9VXVeJzyx3E+pXfv3lZiYuJZ1y0tLbV69epltWnTxvL19bUkWZKs+Ph469ixYx55r4CmxGZZHrxQDaBRS0pK0po1a/TNN980dCmQ9NRTT+mjjz7SrFmzFBAQoNLSUn3xxReaPn26tm7detolRqC54/IQADSQ0aNHa+vWrRozZoyOHj0qHx8fXXzxxXrxxRcJLMAZcKYFAAAYgU8PAQAAIxBaAACAEQgtAADACIQWAABghEb/6aGKigodPHhQAQEBfNcGAACGsCxLv/76q8LDwyt9NcW5aPSh5eDBg4qIiGjoMgAAQC3s379fF1xwgUf21ehDS0BAgKSTLzowMLCBqwEAADVRWFioiIgI13HcExp9aDl1SSgwMJDQAgCAYTw5tYOJuAAAwAiEFgAAYARCCwAAMEKjn9MCACaxLEtlZWUqLy9v6FKAOtWiRQvZ7fZ6vR0JoQUAPKSkpESHDh3SsWPHGroUoF74+fkpLCxM5513Xr30R2gBAA+oqKhQZmamWrRoofDwcJ133nncEBNNlmVZKikpUV5enjIzM9WlSxeP3UCuKoQWAPCAkpISVVRUKCIiQn5+fg1dDlDnHA6HvL29tXfvXpWUlMjX17fO+2QiLgB4UH38bxNoLOr7952/LgAAYAQuDwFAHRuXklZvfb06tme99QXUN860AEAz9+2336pHjx4KCAhQv379dODAgTrvMyUlRX369Knzfvr06SObzSaHw6GePXvqgw8+OG2dsWPHKikpqc5rqaq/Pn36KCUlpU76S0pK0tixY+tk3/WN0AIAzdjx48d14403auLEifrxxx8VEBCgSZMmNXRZHvXXv/5V27Zt07XXXqubbrpJX3zxRaX2hQsX6vHHH3drn2vXrlVkZGSt6qlNf1XZs2dPlZ9Ue/zxx7Vw4UKP9deQCC0A0Iz9+OOPOnr0qMaPH6+IiAglJiY2uY9qOxwOdenSRTNnztStt96qF154oVK7n59fvXzypaH68/X1bTKfaCO0AEAzFhERIZvNplmzZqmsrEzR0dF69913Xe0rV67URRddJH9/f11zzTU6ePCgJCkyMlITJkxQaGioHnvsMd14440KCQnR5s2blZSUpEGDBunqq69WUFCQRowYocLCwhrVs3TpUnXp0kVBQUEaM2ZMpRv1zZ8/X+Hh4QoICNCIESN04sQJt1/vwIEDtXHjxkrPnelyjWVZevTRRxUSEqLg4GBNmjRJlmUpOztbNptNffv21d69e2Wz2WSz2ZSdne3a1maz6YcfflBCQoJatWqlgoKCavuTpM2bN6tz585q27atnn766bOu//tLa76+vurYsaOrX5vNpm+++abSfs92eWjBggWKjIxUeHi4kpKSVFFR4epv2rRpuu+++9SyZUv96U9/0o4dO874ftY3JuIC9cjTEzKZdIlz1bZtW73++usaP368Fi1apJkzZ2r06NGSpCNHjui2227TSy+9pIEDB2rSpEmaNWuW61LDr7/+quTkZCUkJOj999/Xb7/9pjVr1kiS1qxZo6VLl+rll1/W6NGjNW3atNPOcPzRl19+qb/97W9atmyZAgMDdf/99+vZZ5/V9OnTtWPHDj3yyCNau3atQkNDNXLkSC1atMjtS1lhYWHKy8urdr2PP/5YixYt0rp16+Tl5aVBgwbp2muv1fXXX6+jR4/qq6++0r333qtt27ZJkpxOZ6Xtx48fr7i4OL333nvy9/evUW1vvvmmVq1apbKyMg0dOlQxMTEaNGhQldvk5ORo7969uvTSS3X06FFJUkBAQLV9vfPOO5oxY4ZWrFihwMBA3XbbbQoODtaDDz4oSfrHP/6hcePG6fvvv9ddd92lp59+WkuXLq3R66hLhBYAaOZuvvlm9e/fX/Pnz9c999yjbdu26ZlnnlFAQID27t0rp9OpTZs2qaSkpNL/uEePHi1fX1+df/75uvHGG/Xee++ptLRUknTVVVdp1KhRkqTHHntM//u//1ttaFm6dKl27dql/v37Szp5w75TTt0mvri4WJGRkUpLq91/AGw2myzLqnY9Hx8fVVRUqKSkRDExMdq7d6+rLSgoSC1btpSXl5eCgoLOuH337t31zDPPuFXb+PHj9ec//1mSdPvtt+vtt9+uNrQ4nU4FBga66qqpV155RZMnT3adsZkxY4ZmzpzpCi0RERGaPXu2JGnEiBFavny5W6+lrnB5CACasYMHD2rXrl1yOp1KSkrSv//9bz333HPat2+fJOnJJ59Uu3bt9Pjjj5/2RZCn5mWcaX5GRESE6+fw8HDl5ORUW0tWVpbuv/9+bdmyRVu2bNGPP/6o1157TZLUsWNHvfTSS3r44YfVunVr3XrrrcrNzXX79ebk5Oj888+vdr0+ffro0Ucf1ahRoxQSEqK7775bRUVFNe7ngQcecLu2U5d5pJPv36FDh864njt1nM2+ffvUqVMn13KnTp0qBbPevXu7fj7vvPNqFPTqA6EFAJqxt956S+PHj3ct/+Uvf5HdbtfRo0e1bNkyffbZZ9q7d6+++uorDRkypMb7zczMdP28b98+hYWFVbvNBRdcoF9++UWRkZGKjIzUgQMHXJckDh48qJ49eyojI0N79uzR0aNHNXPmTDde6Ukff/yxrrrqqmrX2717t4YPH67t27fru+++07fffquXXnrJ1e7l5VXlgbyml4R+71RQlKQDBw64wpXNZqsUFjdt2lRpu1N3pXUnWERERGjXrl2u5V27dql9+/au5VNnbxobQgsANGPXXHONvv76ay1fvlwHDhxQUlKSQkNDdfHFF7smzx45ckT//ve/NXPmzBofGL/55hu99tpr2rlzp+bOnathw4ZVu82oUaP01ltvadWqVdq1a5emTJniquG7777TgAEDtH79eteZhlMTR6tz/Phx/fzzz0pOTlZqaqoefvjharf55JNPdNNNNykjI0PFxcWn9depUycdPHhQGRkZ+vnnn5WRkVGjWqryyiuv6Ntvv9XXX3+tZcuW6cYbb5QktWvXTl9//bUsy9KGDRsqTZSWTs7T8fPz0+rVq7V3797TJuKeyfjx4/X8889r3bp1ysjIUFJSkhISEs75NdQ15rQAQB1rzBOmu3fvriVLligxMVEHDhxQTEyM3n//ffn4+Gj06NFatWqVunbtqm7duikhIUELFy6s0ad2hgwZosWLF2vSpEkaPHiwEhMTq92md+/eWrBggaZMmaJffvlFN998s6ZNmybp5Kd+EhISdMsttyg/P19XXnmlnnzyyRq9xieeeELJycm6/PLL9fHHH6tHjx7VbnPXXXcpIyNDAwcO1PHjxzVgwABNnDjR1X7BBRdo9uzZ6tevn4qLizVnzhzFxMTUqJ6zGTJkiO68804dPnxYEydO1E033SRJuvfee/XJJ5+oc+fOuvjii/Xkk0/qww8/dG3n7e2tl19+WQkJCTp69KgeeOABXXHFFVX2dcstt+jQoUMaPXq0SkpKdPfdd9fqklZ9s1mN5ULVWRQWFsrpdKqgoKDRnq4CaopPDzVdJ06cUGZmpjp27Fiv9+BojJKSkrRnz546u8MrGo+qfu/r4vjN5SEAAGAELg8BADyqPr/HB80LZ1oAAIARCC0A4EGNfJog4FH1/ftOaAEAD/D29pakSt+VAzR1p37fT/3+1zXmtACAB7Ro0UJBQUGuu7T6+fk1uW9LBk6xLEvHjh1Tbm6ugoKC1KJFi3rpl9ACAB4SGhoqSbW6vTxgoqCgINfvfX0gtACAh9hsNoWFhalt27auLw4Emipvb+96O8NyCqEFADysRYsW9f6POdAcMBEXAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiBO+ICcBmXkubR/b06tqdH9wegeeNMCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAj8C3PgME8/a3MANCYcaYFAAAYgdACAACMQGgBAABGcDu0fPfdd/rzn/+sgIAADRw4UPv27ZMkpaenKzo6Wg6HQwMGDFBubq5rm6raAAAAasLt0HLTTTdp8ODB2rFjhyIjI3XXXXepoqJCw4cP1+DBg7Vz5045HA5NnjxZkqpsAwAAqCmbZVlWTVfOy8tT27ZtdejQIYWGhmrDhg3q16+fVq9eraFDh+rw4cOy2+1KT09XXFyc8vLytHHjxrO2+fv7V9tnYWGhnE6nCgoKFBgYeE4vFmhoze3TPq+O7dnQJQBoIHVx/HbrTEtwcLAuuOACffTRR5KkNWvWKDo6WuvXr1dsbKzs9pOfoI6OjlZ5ebnS09OrbAMAAKgpt+7TYrfbtWLFCvXt21f33HOPWrZsqY0bN2r+/Plq06aNaz0vLy8FBwcrJydH2dnZZ207k+LiYhUXF7uWCwsL3X1NAACgCXLrTMvx48c1ZswYzZgxQ+np6Ro1apTuuusuSdIfrzJZliWbzVZt2x/Nnj1bTqfT9YiIiHCnRAAA0ES5dabl448/VklJiR577DFJ0jPPPCN/f39dffXV2r59u2u98vJy5efnKzQ0VGFhYWdtO5OpU6dqypQpruXCwkKCCxpMc5uDAgCNmVtnWlq0aKHy8nLXsmVZqqioUN++fZWWlqaysjJJUkZGhux2u2JiYhQfH3/WtjPx8fFRYGBgpQcAAIBboaVXr14qLCzU/PnzlZWVpSeffFIRERHq1auXQkJClJiYqKysLCUnJ2vYsGHy8/NTXFzcWdsAAABqyq3QEhISotTUVC1evFgXXXSRvvrqK7377rvy8fFRamqqVq9eraioKJ04cULz5s072YGX11nbAAAAasqt+7Q0BO7TgobEnJZzw31agOarwe/TAgAA0FAILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBHsDV0A4EnjUtIaugQAQB3hTAsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEewNXQCApmtcSppH9/fq2J4e3R8As3CmBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxgb+gCAKCmxqWkeXR/r47t6dH9AahbnGkBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEt0NLaWmpEhISFBAQoEsuuUQbN26UJKWnpys6OloOh0MDBgxQbm6ua5uq2gAAAGrC7dDy7LPPas+ePcrIyNCIESM0atQoVVRUaPjw4Ro8eLB27twph8OhyZMnS1KVbQAAADVlsyzLcmeDqKgovfPOO7r00kv122+/ac2aNQoODtawYcN0+PBh2e12paenKy4uTnl5edq4caOGDh16xjZ/f/9q+yssLJTT6VRBQYECAwNr/ULRPHj6Ph5o2rhPC1B36uL47daZluzsbO3evVvr1q2T0+lU7969demll2rDhg2KjY2V3X7yXnXR0dEqLy9Xenq61q9ff9Y2AACAmnIrtBw6dEheXl769ttvtXXrVnXt2lUTJkxQdna22rRp89+denkpODhYOTk5VbadSXFxsQoLCys9AAAA3AotRUVFKi8vV2JioiIjIzVp0iR9/vnnqqio0B+vMlmWJZvN5vr5bG1/NHv2bDmdTtcjIiLCnRIBAEAT5VZocTqdkqRWrVpJklq3bi3LstSuXTvl5eW51isvL1d+fr5CQ0MVFhZ21rYzmTp1qgoKClyP/fv3u/2iAABA0+NWaImKipK3t7f+85//SJJycnLUokULxcfHKy0tTWVlZZKkjIwM2e12xcTEVNl2Jj4+PgoMDKz0AAAAcCu0OBwO3XDDDUpKStKuXbv0wgsv6LrrrlNcXJxCQkKUmJiorKwsJScna9iwYfLz86uyDQAAoKbcvk/LwoULZVmWunXrppycHP3973+Xl5eXUlNTtXr1akVFRenEiROaN2/eyQ6qaAMAAKgpt+/TUt+4TwvcwX1a4A7u0wLUnQa/TwsAAEBDIbQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEawN3QBaN64gy0AoKY40wIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBHsDV0AzDIuJa2hSwA8xtO/z6+O7enR/QGojDMtAADACIQWAABgBEILAAAwQq1CyxdffCGbzaa1a9dKktLT0xUdHS2Hw6EBAwYoNzfXtW5VbQAAADXldmgpLS3Vvffe61quqKjQ8OHDNXjwYO3cuVMOh0OTJ0+utg0AAMAdbn96aP78+Wrbtq2ysrIkSevWrdORI0eUlJQku92uxMRExcXFqaioSBs3bjxrm7+/v8dfDAAAaLrcOtOSlZWlOXPmaMGCBa7n1q9fr9jYWNntJ/NPdHS0ysvLlZ6eXmXb2RQXF6uwsLDSAwAAwK3Q8uCDDyohIUFdu3Z1PZedna02bdr8d4deXgoODlZOTk6VbWcze/ZsOZ1O1yMiIsKdEgEAQBNV49CyZs0abd68WdOmTTutzbKs05ZtNlu1bWcydepUFRQUuB779++vaYkAAKAJq3Foeeutt3To0CGFh4crKChIBQUFGjx4sMLCwpSXl+dar7y8XPn5+QoNDa2y7Wx8fHwUGBhY6QEAAFDj0PLMM89ox44d2rJli7Zs2aKAgAAtWrRI8fHxSktLU1lZmSQpIyNDdrtdMTExVbYBAAC4o8ahpU2bNoqMjHQ9vLy8FBoaqri4OIWEhCgxMVFZWVlKTk7WsGHD5OfnV2UbAACAO875jrheXl5KTU3V6tWrFRUVpRMnTmjevHnVtgEAALij1t/ynJ+f7/r5sssu09atW8+4XlVtAAAANcV3DwEAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMIK9oQsAgKZiXEqaR/f36tieHt0fYDrOtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGsDd0Aahb41LSGroEAAA8gjMtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYAS3Q8vu3bvVu3dvBQQEqE+fPtq7d68kKT09XdHR0XI4HBowYIByc3Nd21TVBgAAUBM2y7Isdzbo16+fwsLCNGvWLE2ZMkXFxcVatWqVOnfurDvuuEMTJkzQfffdp5YtW+qNN95QRUXFWdtqorCwUE6nUwUFBQoMDKzVi2zOxqWkNXQJABqJV8f2bOgS0IzUxfHbrdBSUlIiX19fff/997rkkkv04YcfauTIkXr//fc1dOhQHT58WHa7Xenp6YqLi1NeXp42btx41jZ/f/9q+yS0nBtCC4BTCC2oT3Vx/Hbr8lBpaanmzp2rjh07SpIOHz4sh8Oh9evXKzY2Vna7XZIUHR2t8vJypaenV9l2JsXFxSosLKz0AAAAcCu0+Pv76+GHH5bD4VBpaalefPFFjR49WtnZ2WrTps1/d+rlpeDgYOXk5FTZdiazZ8+W0+l0PSIiImr50gAAQFNSq08PlZWV6Y477pCXl5eSk5MlSX+8ymRZlmw2W7VtfzR16lQVFBS4Hvv3769NiQAAoImxu7tBRUWFRowYod27d+vTTz+Vw+FQWFiYtm/f7lqnvLxc+fn5Cg0NrbLtTHx8fOTj41OLlwIAAJoyt8+0JCcn6+eff9Znn32mVq1aSZLi4+OVlpamsrIySVJGRobsdrtiYmKqbAMAAKgpt0JLdna25s+fr//7v/+TJOXn5ys/P19xcXEKCQlRYmKisrKylJycrGHDhsnPz6/KNgAAgJpyK7R89NFHKiws1JVXXqng4GDXY9++fUpNTdXq1asVFRWlEydOaN68eSc78PI6axsAAEBNuX1zufrGfVrODfdpAXAK92lBfWrw+7QAAAA0FEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYAR7QxeAysalpDV0CQAANEqcaQEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEfjuIQBoJjz53Wavju3psX0BNcWZFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEbhPCwDAbZ6854vEfV9QM5xpAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBHsDV2A6calpDV0CQAANAucaQEAAEYgtAAAACMQWgAAgBGY0wIAaHCenh/46tieHt0fGgfOtAAAACMQWgAAgBEILQAAwAj1ElrS09MVHR0th8OhAQMGKDc3tz66BQAATYjNsiyrLjuoqKhQ586ddccdd2jChAm677771LJlS73xxhs12r6wsFBOp1MFBQUKDAw853q4GRwAwF1M7HWfp4/fUj18emjdunU6cuSIkpKSZLfblZiYqLi4OBUVFcnf37+uuwcAAE1EnYeW9evXKzY2Vnb7ya6io6NVXl6u9PR0xcfHn7Z+cXGxiouLXcsFBQWSTiY2Tyg5/ptH9gMAaD48dQxqTk69Z568oFPnoSU7O1tt2rRxLXt5eSk4OFg5OTlnXH/27NmaMWPGac9HRETUWY0AAFTl9XsbugJzHT58WE6n0yP7qpeby/0xZVmWJZvNdsZ1p06dqilTpriW8/Pz1aFDB+3bt89jLxq1U1hYqIiICO3fv99j1ydRO4xF48FYNC6MR+NRUFCg9u3bq1WrVh7bZ52HlrCwMG3fvt21XF5ervz8fIWGhp5xfR8fH/n4+Jz2vNPp5BewkQgMDGQsGgnGovFgLBoXxqPx8PLy3AeV6/wjz/Hx8UpLS1NZWZkkKSMjQ3a7XTExMXXdNQAAaELqPLTExcUpJCREiYmJysrKUnJysoYNGyY/P7+67hoAADQhdR5avLy8lJqaqtWrVysqKkonTpzQvHnzary9j4+PEhMTz3jJCPWLsWg8GIvGg7FoXBiPxqMuxqLOby4HAADgCXz3EAAAMAKhBQAAGIHQAgAAjEBoAQAARmjw0JKenq7o6Gg5HA4NGDBAubm51W7z0Ucf6cILL5S/v79GjBihY8eO1UOlTV9txmLz5s267LLLFBgYqBtvvFFHjhyph0qbvtqMxSlLly6VzWbTnj176q7AZqQ2Y1FYWKjhw4fL399fvXr10s8//1wPlTZ9tRmLdevWqXv37goMDNStt96qo0eP1kOlzcehQ4fUu3dvbdmypUbrn+vxu0FDS0VFhYYPH67Bgwdr586dcjgcmjx5cpXb5Ofn69Zbb9WUKVP0ww8/aM+ePXr66afrqeKmqzZjUVFRodtvv10DBgzQtm3blJubq+nTp9dTxU1XbcbilPz8fD3yyCN1XGHzUduxeOSRRxQQEKAffvhBl156qe69ly+uOVe1GYtff/1VQ4cO1UMPPaQffvhBJ06c4O/DgxISEhQeHq4vvviiRut75PhtNaDPPvvMCgwMtEpLSy3LsqzNmzdbDofD+u233866zeLFi61LLrnEtfzOO+9Y7du3r/Nam7rajMXOnTstSa51FixYYHXv3r1e6m3KajMWp0ycONEaOnSoJcnKzMys40qbvtqMxbFjx6zg4GDrl19+sSzLsrKzs60PPvigXuptymozFhs3brR8fHxcy8uWLat0/MC5ycvLszIzMy1JVkZGRrXre+L43aBnWtavX6/Y2FjZ7Se/Aik6Olrl5eVKT0+vcpsrr7zStdyrVy/t27dP+/fvr/N6m7LajIWfn5+ef/55+fv7Szr5TZ4Oh6Ne6m3KajMW0slLdampqZo7d259lNks1GYstm7dqoCAAD333HNq2bKlRo4cqdjY2PoqucmqzVh06tRJdrtd69atU3l5uT755BNFR0fXU8VNX5s2bRQZGVnj9T1x/G7Q0JKdna02bdq4lr28vBQcHKycnJwab9O6dWtJqnIbVK82YxEeHq4HH3xQ0slv81y0aJFGjx5d57U2dbUZi4qKCk2cOFGzZs1SSEhIfZTZLNRmLA4dOqS8vDwVFRXpu+++k81m09SpU+uj3CatNmPRunVrpaSkqG/fvvLz89PatWv14osv1ke5OANPHL8bfCKu9Ycb8lqWJZvNVuNtTv1c3TaoXm3GQpKKioo0ZMgQxcTEaOLEiXVVXrPi7lj84x//kM1m0913313XpTU77o5FUVGRSkpKNGfOHHXs2FF333231qxZU9dlNgvujsWhQ4c0adIkLV68WJs2bVLPnj1rPD8MdeNcj98NGlrCwsKUl5fnWi4vL1d+fr5CQ0NrvM3hw4clqcptUL3ajIUkHTt2TAMHDpSfn5/efPNNj34FeXNVm7FYsWKFtm3bplatWqlDhw6SpO7du+urr76q83qbstqMhdPplK+vr+tSaevWrV3/TqH2avt3cdFFF2ns2LHq1q2b5syZo9dff10FBQX1UTL+wBPH7wY9wsTHxystLU1lZWWSpIyMDNntdsXExFS5zddff+1a/uabbxQZGal27drVeb1NWW3GQpLuvfde+fv7a+XKlfL19a2PUpu82ozF8uXLtX37dm3ZskVffvmlJOnDDz/U5ZdfXi81N1W1GYuuXbuqqKhIBw8elHTy1Df/qTp3tRmLFi1aqLy83LVcVlYmy7Jc82JQvzxy/HZr2q6HlZeXW506dbKeeOIJa//+/daQIUOsUaNGWZZlWQUFBVZJSclp2+Tn51tOp9NauHChlZmZacXGxlpPPfVUfZfe5NRmLLZs2WI5nU7rP//5j3X06FHXA+emNmPxe0ePHuXTQx5S27Ho0aOHNWbMGGvXrl3WVVddZU2aNKk+y26SajMW33//veXt7W0tX77c2rt3rzVy5EjryiuvrO/Smzz94dNDdXn8btDQYlknP7bWvXt3y8fHx+rfv7+Vm5trWZZldejQwXrvvffOuM1HH31kRUVFWQ6Hw7rtttusoqKieqy46XJ3LJKSkixJpz1w7mrzd3EKocWzajMWO3bssHr27Gn5+/tbgwcPto4cOVKPFTddtRmLf/7zn1ZUVJTl7+9v9evXz9q1a1c9Vtw8/DG01OXx2/b/OwQAAGjUmDUJAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEb4f0iLsLWE8vUGAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.stats import beta, binom\n",
    "\n",
    "\n",
    "class MetropolisHastings:\n",
    "    def __init__(self, proposal_dist, accepted_dist, m=1e4, n=1e5):\n",
    "        \"\"\"\n",
    "        Metropolis Hastings\n",
    "\n",
    "        :param proposal_dist: 建议分布\n",
    "        :param accepted_dist: 接受分布\n",
    "        :param m: 收敛步数\n",
    "        :param n: 迭代步数\n",
    "        \"\"\"\n",
    "        self.proposal_dist = proposal_dist\n",
    "        self.accepted_dist = accepted_dist\n",
    "        self.m = m\n",
    "        self.n = n\n",
    "\n",
    "    @staticmethod\n",
    "    def __calc_acceptance_ratio(q, p, x, x_prime):\n",
    "        \"\"\"\n",
    "        计算接受概率\n",
    "\n",
    "        :param q: 建议分布\n",
    "        :param p: 接受分布\n",
    "        :param x: 上一状态\n",
    "        :param x_prime: 候选状态\n",
    "        \"\"\"\n",
    "        prob_1 = p.prob(x_prime) * q.joint_prob(x_prime, x)\n",
    "        prob_2 = p.prob(x) * q.joint_prob(x, x_prime)\n",
    "        alpha = np.minimum(1., prob_1 / prob_2)\n",
    "        return alpha\n",
    "\n",
    "    def solve(self):\n",
    "        \"\"\"\n",
    "        Metropolis Hastings 算法求解\n",
    "        \"\"\"\n",
    "        all_samples = np.zeros(self.n)\n",
    "        # (1) 任意选择一个初始值\n",
    "        x_0 = np.random.random()\n",
    "        for i in range(int(self.n)):\n",
    "            x = x_0 if i == 0 else all_samples[i - 1]\n",
    "            # (2.a) 从建议分布中抽样选取\n",
    "            x_prime = self.proposal_dist.sample()\n",
    "            # (2.b) 计算接受概率\n",
    "            alpha = self.__calc_acceptance_ratio(self.proposal_dist, self.accepted_dist, x, x_prime)\n",
    "            # (2.c) 从区间 (0,1) 中按均匀分布随机抽取一个数 u\n",
    "            u = np.random.uniform(0, 1)\n",
    "            # 根据 u <= alpha，选择 x 或 x_prime 进行赋值\n",
    "            if u <= alpha:\n",
    "                all_samples[i] = x_prime\n",
    "            else:\n",
    "                all_samples[i] = x\n",
    "\n",
    "        # (3) 随机样本集合\n",
    "        samples = all_samples[self.m:]\n",
    "        # 函数样本均值\n",
    "        dist_mean = samples.mean()\n",
    "        # 函数样本方差\n",
    "        dist_var = samples.var()\n",
    "        return samples[self.m:], dist_mean, dist_var\n",
    "\n",
    "    @staticmethod\n",
    "    def visualize(samples, bins=50):\n",
    "        \"\"\"\n",
    "        可视化展示\n",
    "        :param samples: 抽取的随机样本集合\n",
    "        :param bins: 频率直方图的分组个数\n",
    "        \"\"\"\n",
    "        fig, ax = plt.subplots()\n",
    "        ax.set_title('Metropolis Hastings')\n",
    "        ax.hist(samples, bins, alpha=0.7, label='Samples Distribution')\n",
    "        ax.set_xlim(0, 1)\n",
    "        ax.legend()\n",
    "        plt.show()\n",
    "\n",
    "        \n",
    "class ProposalDistribution:\n",
    "    \"\"\"\n",
    "    建议分布\n",
    "    \"\"\"\n",
    "\n",
    "    @staticmethod\n",
    "    def sample():\n",
    "        \"\"\"\n",
    "        从建议分布中抽取一个样本\n",
    "        \"\"\"\n",
    "        # B(1,1)\n",
    "        return beta.rvs(1, 1, size=1)\n",
    "\n",
    "    @staticmethod\n",
    "    def prob(x):\n",
    "        \"\"\"\n",
    "        P(X = x) 的概率\n",
    "        \"\"\"\n",
    "        return beta.pdf(x, 1, 1)\n",
    "\n",
    "    def joint_prob(self, x_1, x_2):\n",
    "        \"\"\"\n",
    "        P(X = x_1, Y = x_2) 的联合概率\n",
    "        \"\"\"\n",
    "        return self.prob(x_1) * self.prob(x_2)\n",
    "\n",
    "class AcceptedDistribution:\n",
    "    \"\"\"\n",
    "    接受分布\n",
    "    \"\"\"\n",
    "\n",
    "    @staticmethod\n",
    "    def prob(x):\n",
    "        \"\"\"\n",
    "        P(X = x) 的概率\n",
    "        \"\"\"\n",
    "        # Bin(4, 10)\n",
    "        return binom.pmf(4, 10, x)\n",
    "\n",
    "    \n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# 收敛步数\n",
    "m = 1000\n",
    "# 迭代步数\n",
    "n = 10000\n",
    "\n",
    "# 建议分布\n",
    "proposal_dist = ProposalDistribution()\n",
    "# 接受分布\n",
    "accepted_dist = AcceptedDistribution()\n",
    "\n",
    "metropolis_hastings = MetropolisHastings(proposal_dist, accepted_dist, m, n)\n",
    "\n",
    "# 使用 Metropolis-Hastings 算法进行求解\n",
    "samples, dist_mean, dist_var = metropolis_hastings.solve()\n",
    "print(\"均值:\", dist_mean)\n",
    "print(\"方差:\", dist_var)\n",
    "\n",
    "# 对结果进行可视化\n",
    "metropolis_hastings.visualize(samples, bins=20)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d67c0b6d",
   "metadata": {},
   "source": [
    "## 习题19.8\n",
    "![image.png](./images/exercise8.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a5cdaed2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theta均值：0.5213045386327443, theta方差：0.018263428918901512\n",
      "eta均值：0.12231449922210416, eta方差：0.006411499241635615\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAAGvCAYAAABMwk8eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAtMklEQVR4nO3de3TU9Z3/8dcMA2ESkhBIaIIGiEYRFmGCEFsMjbIathYtwraAiodKFW+7HrRuRVcTozZeKlirtFK11Aos5ai1phbWVUSJrolOgoKIQRGSlVwAk5EAucx8fn/wY2wgCblMMvmE5+Oc7zn5fj/fy3vyMczLz/fmMMYYAQAAWMwZ7gIAAAC6ikADAACsR6ABAADWI9AAAADrEWgAAID1CDQAAMB6BBoAAGA9Ag0AALAegQYAAFiPQANY6plnnlFKSorcbrd++MMf6quvvjphnbfeeksOh+OE5Tk5Ofrud7/b6r6//PJLORwOffrppyGtuSW/+c1vlJqaqoEDByojI0Mff/xxtx+zNQ6HQ2+99Var7RdeeKFycnJ6rB4A7UegASz08ssv67rrrtOCBQv04osvqry8XFdfffUJ602cOFHvvfdeGCpsn6efflp33XWX7rzzTr388ss6dOiQZs+eraampnCX1qLly5frZz/7WbjLANACV7gLANBxOTk5ysrKUnZ2tiTJ7XZr2rRp+uijjzR+/PjgejExMW2OxITb008/rQULFgRDQnx8vNLT01VQUKDMzMwwV3eisWPHhrsEAK1ghAawzFdffaWPPvpIWVlZwWUTJ06UJG3dujVcZXXK/v37deDAgeD8hAkT9Prrr2v06NFhrAqAjQg0gGW2b98uSUpJSQkui42N1fbt2/WDH/yg2bqtXUNzzD333KPBgwdr2LBhys3NlTGmWfs777yjc845J3idzv/93/81a3/xxRd17rnnyu12KyUlRU8++WSHPstll12mNWvW6Pbbb9e+ffs0YMAAXXzxxUpMTAyus2rVKo0dO1aRkZE655xztGbNmmCbw+FQXl6eRo4cqdNPP12bNm3SuHHjNGTIEL366qvBdX75y19q/PjxGjRokP75n/9ZO3bs6FCdx7R2Dc2oUaO0cuVK5eXlKTExUXFxcfq3f/u3Zr/Pjz/+WFOmTFF0dLQuueQS5eXlKTk5WX/60586VQuA4xgAVlmzZo2RZF5//XVjjDFNTU2msbHRNDY2mkAg0GzdjRs3mpb+zLOzs43L5TLf+973zN/+9jdz//33G6fTaVasWGGMMWbXrl1GkomLizMrVqwwL774oklJSTFTpkwJ7qO0tNT069fPXH/99ebNN980eXl5RpIpKCho92c5dOiQmT9/vpFkIiMjzd13320OHz4cbN+8ebNxOBzmtttuMxs3bjT33nuvcblc5vPPPzfGGCPJnHPOOeZvf/ubSUhIMIMGDTKrVq0y06ZNM5deemlwnYEDB5rHHnvM5Ofnm/T0dJOcnGwOHTp0Qj2SzMaNG1utNzMz02RnZ5+wfOTIkcbj8ZhJkyaZv/zlLyY3N9dIMn/961+D66Smppr58+ebN954w0yfPt2cffbZ5t133zW7d+9u9+8LQOsINIBlnn/+eSPJvPHGG8YYY0aPHm0kGUnmD3/4Q7N12wo0/fv3N+Xl5cFlc+fONWPGjDHGfBtoHn744WD7+vXrjSTz/vvvN9v3sXljjHnttdfMnj17OvyZPvjgA/PDH/7QSDIXXHCBOXLkiDHGmA8//NCsWLHCNDU1GWOMqa6uNv379zdr1641xhwNIC+88IIxxpjvf//75mc/+5kxxph7773XZGZmBtf593//9+CxSktLjSSzatWqE+roSqBJTEw0Pp8vuGzs2LHmgQceCNYtyXzyySfGmKO/pwEDBrTnVwOgnTjlBFhm0KBBkqSDBw9KOnrap6ioSElJSR3az5lnnqnTTjstOP/d735XO3fuVCAQCC77xwtzj11cfOx0zZQpUzR58mRdeumlmj9/vp588klNmDBBycnJHf5M5513nvLz87V8+XIVFBToqaeeknT02qBzzz1Xd955p6ZMmaLTTz9dTU1NOnToUHDbY5/B4XA0+/kfXXDBBcGfU1NTNWTIEJWWlna4zrb89Kc/VXR0dHA+ISFBjY2NkqShQ4cqISFBr7zyinw+n/Lz8zVmzJiQHh841RFoAMuceeaZkqRdu3ZJkv7pn/5JkyZNUn19fYf243Q2//Pv16+fjDHNAs0/rtOvXz9Jkt/vlyQNGDBA7777rtasWaOUlBQ9++yzOuuss/TBBx+06/her1fjxo3Ttm3bgstuvPFGjRs3Tl6vV9LRu6C+//3vq6amRjfccIM++eQTjRgxokOfU9IJ1wYFAoETPn9XHeuX1o6flpam+++/X7GxsXrppZf0u9/9LqTHB051BBrAMuPGjVNiYqL+8pe/BJft2bOn2d1C7bFz505VVlYG599//32NGjVKLte3T3PYvHlzs3ZJOuussyRJGzdu1GOPPaZLLrlEubm5+uCDDxQXF9fui1wTExO1bdu24H4lqampSV9//bVGjhwpSfrd736nuXPn6ve//72uueYaDRo0qMOfU5I2bdoU/PnTTz9VTU2Nzj777A7vpy3HAl9LXnnlFR04cEDV1dXavn27du/e3atvpwdsxHNoAMs4nU7l5ubq+uuv16JFizRz5kw9+uijHd5PU1OTfvzjH+uuu+5SSUmJ1qxZo2XLljVb59iIQlxcnP7jP/5D5513nqZMmSLp6AjNPffcI5fLpcmTJ2vbtm2qrKzUGWec0a7jDx8+XD/60Y/0i1/8QsYYjRgxQn/4wx/09ddfa+HChZKOnqp57733tGHDBlVUVOjBBx/UN9980+EH761cuVJnnHGGxowZo+zsbI0cOVIzZ87s0D66wu1269NPP9Xq1as1evRoHT58WKeffroSEhJ6rAagryPQABa67rrr1NDQoF/+8pd66aWXNHPmzA6PXEyaNElpaWn68Y9/rLi4ON133326+eabm61z77336oEHHtC+fft00UUXafny5cHrUy644AKtWLFCjz32mO655x7FxsbqlltuOWEfbXn++ed15513KicnRz6fTx6PR2+88UYwFP3mN7/R9ddfr5kzZyoxMVG33nqrVq1apXfeeadDT+y9//779cILL+iTTz7RpEmTlJ+fr4iIiHZv31UZGRkaOnSo7r77btXW1gZPD/7rv/6r1q1b12N1AH2Zwxx/chkA+hCHw6G///3v+pd/+Zew1XD11VfrwIEDuu222xQVFaXDhw/r5Zdf1vLly7V//34NHjw4bLUBfQUjNADQzW655RbdfffdmjNnjnw+nyIjIzV+/Hi98MILhBkgRBihAQAA1uMuJwAAYD0CDQAAsB6BBgAAWI9AAwAArGf1XU6BQEBfffWVoqOjT3h3CwAA6J2MMfrmm280fPjwkL2GxOpA89VXX3XqRXgAACD8ysrKdPrpp4dkX1YHmmNvti0rK1NMTEyYqwEAAO3h8/mUnJzc7A31XWV1oDl2mikmJoZAAwCAZUJ5uQgXBQMAAOsRaAAAgPUINAAAwHpWX0MDAIBNjDFqamqS3+8Pdyndql+/fnK5XD36SBUCDQAAPaChoUF79+7VoUOHwl1Kj4iMjFRSUpIGDBjQI8cj0AAA0M0CgYB27dqlfv36afjw4RowYECffSCsMUYNDQ2qrq7Wrl27dNZZZ4Xs4XltIdAAANDNGhoaFAgElJycrMjIyHCX0+3cbrf69++v3bt3q6GhQQMHDuz2Y3JRMAAAPaQnRip6i57+rKfObxYAAPRZBBoAAGC9DgeavXv3KjMzUyUlJSe0vf3223I4HHrrrbeCy7xerzwej9xut7KyslRVVdWuNgAAgPbqUKBZtGiRhg8frrfffvuEtsbGRt10003NlgUCAc2ePVszZsxQaWmp3G63Fi9efNI2AACAjujQXU4PPviglixZopSUlBPali1bpmHDhqm8vDy4bNOmTTpw4IBycnLkcrmUnZ2tjIwM1dXVqbCwsNW2qKiorn8yAAAssHBlUY8e79kFkzu97a9//Ws9/vjjqqysVFZWlv74xz8qNjY2hNV1XodGaOLj4zVq1KgTlpeXl+uhhx7SU0891Wx5QUGB0tPT5XIdzU0ej0d+v19er7fNttbU19fL5/M1mwAAQPe766679OSTT+qPf/yjNm/erOLiYt13333hLisoJM+hufXWW7Vo0SKNGTOm2fKKigrFx8cH551Op+Li4lRZWdlmW2vy8vJ6zy9v9ZzOb3vl2tDVAQBANysqKtLDDz+soqIiTZw4UdLRy1BWrlyppUuXhrm6o7p8l9P69ev14Ycf6p577mmx3RhzwvyxpyO21daSJUuWqLa2NjiVlZV1sXoAAHAyv/rVrzRt2rRgmJGkhIQE7du3L4xVNdflQLN27Vrt3btXw4cP1+DBg1VbW6sZM2Zo9erVSkpKUnV1dXBdv9+vmpoaJSYmttnWmoiICMXExDSbAABA96mvr9err76qK664otnyw4cP95rrZ6QQBJpHH31UO3bsUElJiUpKShQdHa1nnnlGl19+uaZOnaqioiI1NTVJkoqLi+VyuZSWltZmGwAA6B28Xq8OHz6s22+/XYMGDQpOd9xxh0aPHh3u8oK6fA1NfHz8CdfCJCYmatCgQcrIyFBCQoKys7N14403Kjc3V7NmzVJkZGSbbQAAoHf47LPPNHDgQH388cfNll9++eW64IILwlTVibr1ScFOp1Pr1q1Tfn6+UlNTdeTIkeDFQ221AQCA3sHn82nYsGFKTU0NTgMGDNCnn36q2bNnh7u8oE6N0Bx/Me8/qqmpaTY/ceJEbdmypcV122oDAADhFx8fL5/P1+zGnQcffFCXXnqpxo4dG+bqvhWS27YBAEDfNG3aNB05ckQPPfSQ5s2bp9WrV+uvf/2rCgsLw11aMwQaAADCqCtP7u0J3/nOd7Ry5Urdcccduv/++zVt2jRt3rxZycnJ4S6tGQINAABo05w5czRnThceKNsDuvWiYAAAgJ5AoAEAANYj0AAAAOsRaAAAgPUINAAAwHoEGgAAYD0CDQAAsB7PoQEAIJxW9/DzXa5c27PH6yGM0AAAAOsRaAAAgPUINAAA4KTGjRunBx54QDfccIOGDBmixMRELVu2LNxlBRFoAABAm+rr67Vjxw49//zzyszMVGFhoa666irdeeedOnjwYLjLk0SgAQAAJ7F161Y1NTXpiSee0Lx585Samqprr71WDQ0NOnToULjLk0SgAQAAJ1FSUqLExERNnz49uKyyslIDBgzQ0KFDw1jZtwg0AACgTVu2bNGkSZPkcDiaLRs3bpz69esXxsq+RaABAABtKikpkcfjOemycCLQAACANn300Ue9PtDwpGAAAMKplz+598svv1RtbW2z8NLQ0KDt27crLS0tfIUdh0ADAABaNWrUKBljmi3btm2bmpqaNGHChDBVdSJOOQEAgA4pKSnRGWecoejo6HCXEkSgAQAAHbJly5Zedf2MxCknAADQQY8//ni4SzgBIzQAAMB6BBoAAGA9Ag0AALDeqX0Nzeo54a4AAHAKOf72576spz8rIzQAAHSz/v37S1KveTN1Tzj2WY999u52ao/QAADQA/r166fBgwerqqpKkhQZGdnsRY99iTFGhw4dUlVVlQYPHtxjL68k0AAA0AMSExMlKRhq+rrBgwcHP3NPINAAANADHA6HkpKSNGzYMDU2Noa7nG7Vv3//HhuZOYZAAwBAD+rXr1+Pf9mfCrgoGAAAWI9AAwAArEegAQAA1iPQAAAA63FRcDh05QnFV64NXR0AAPQRHR6h2bt3rzIzM1VSUhJc9sUXXygzM1PR0dG68MILtXv37mCb1+uVx+OR2+1WVlZWs/vv22oDAABorw4FmkWLFmn48OF6++23my2//vrrNWLECG3dulVDhw7VzTffLEkKBAKaPXu2ZsyYodLSUrndbi1evPikbQAAAB3hMB14e9S+fft08OBBpaSkqLi4WB6PRw0NDRo4cKC2bt2qsWPH6rXXXtO8efNUW1urjRs3aubMmdq/f79cLpe8Xq8yMjJUXV2twsLCVtuioqLaVY/P51NsbKxqa2sVExPT8U9v48spOeUEALBcl7+/W9ChEZr4+HiNGjWq2bLGxkY98sgjSklJkSTt379fbrdbklRQUKD09HS5XEcv1fF4PPL7/fJ6vW22AQAAdESXLwqOiorSz3/+c0lHw80TTzyh+fPnS5IqKioUHx8fXNfpdCouLk6VlZVttrWmvr5e9fX1wXmfz9fV8gEAQB8Qstu2m5qadNVVV8npdCo3Nze4/PgzWsaY4BtG22prSV5enmJjY4NTcnJyqMoHAAAWC0mgCQQCmjt3rnbu3Km///3vwVNOSUlJqq6uDq7n9/tVU1OjxMTENttas2TJEtXW1gansrKyUJQPAAAsF5JAk5ubq507d+rNN9/UkCFDgsunTp2qoqIiNTU1SZKKi4vlcrmUlpbWZltrIiIiFBMT02wCAADocqCpqKjQsmXL9Nvf/laSVFNTo5qaGgUCAWVkZCghIUHZ2dkqLy9Xbm6uZs2apcjIyDbbAAAAOqLLgWbDhg3y+XyaMmWK4uLigtOePXvkdDq1bt065efnKzU1VUeOHNHSpUuPHriNNgAAgI7o0HNoehueQwMAgH3C/hwaAACA3ohAAwAArEegAQAA1iPQAAAA6xFoAACA9Qg0AADAegQaAABgPQINAACwHoEGAABYj0ADAACsR6ABAADWI9AAAADrEWgAAID1CDQAAMB6BBoAAGA9Ag0AALAegQYAAFiPQAMAAKxHoAEAANYj0AAAAOsRaAAAgPUINAAAwHoEGgAAYD0CDQAAsB6BBgAAWI9AAwAArEegAQAA1iPQAAAA6xFoAACA9Qg0AADAegQaAABgPQINAACwHoEGAABYj0ADAACsR6ABAADWI9AAAADrEWgAAID1CDQAAMB6BBoAAGA9Ag0AALBehwPN3r17lZmZqZKSkuAyr9crj8cjt9utrKwsVVVVdbkNAACgvRzGGNPelRctWqQVK1ZIkoqLi+XxeBQIBHTmmWfqqquu0g033KCbb75ZgwYN0qpVqzrd1l4+n0+xsbGqra1VTExMxz/96jkd38ZmV64NdwUAAHT9+7sFHQo0+/bt08GDB5WSkhIMNBs3btTMmTO1f/9+uVwueb1eZWRkqLq6WoWFhZ1qi4qKalc9BJoOItAAAHqB7gg0HTrlFB8fr1GjRjVbVlBQoPT0dLlcLkmSx+OR3++X1+vtdFtr6uvr5fP5mk0AAABdvii4oqJC8fHx3+7Q6VRcXJwqKys73daavLw8xcbGBqfk5OSulg8AAPqAkNzldPxZK2OMHA5Hl9pasmTJEtXW1gansrKyUJQPAAAs1+VAk5SUpOrq6uC83+9XTU2NEhMTO93WmoiICMXExDSbAAAAuhxopk6dqqKiIjU1NUk6eveTy+VSWlpap9sAAAA6osuBJiMjQwkJCcrOzlZ5eblyc3M1a9YsRUZGdroNAACgI7ocaJxOp9atW6f8/HylpqbqyJEjWrp0aZfaAAAAOqJDz6HpbXgOTQfxHBoAQC/QHc+hcYVkLwDQAxauLAr5Pp9dMDnk+wTQ83g5JQAAsB6BBgAAWI9AAwAArEegAQAA1iPQAAAA6xFoAACA9Qg0AADAegQaAABgPQINAACwHoEGAABYj0ADAACsR6ABAADWI9AAAADrEWgAAID1CDQAAMB6BBoAAGA9Ag0AALCeK9wFAOi7Fq4sCncJAE4RjNAAAADrMUID4JTWHaNIzy6YHPJ9AmgbIzQAAMB6BBoAAGA9Ag0AALAegQYAAFiPQAMAAKxHoAEAANYj0AAAAOsRaAAAgPUINAAAwHoEGgAAYD0CDQAAsB6BBgAAWI+XUwJAiIX6hZe87BI4OUZoAACA9Qg0AADAegQaAABgPQINAACwHoEGAABYL6SB5uOPP9b3vvc9RUdHa/r06dqzZ48kyev1yuPxyO12KysrS1VVVcFt2moDAABoj5AGmiuuuEIzZszQjh07NGrUKF177bUKBAKaPXu2ZsyYodLSUrndbi1evFiS2mwDAABor5A9h6a6ulqff/65Fi5cqMTERC1YsEAXX3yxNm3apAMHDignJ0cul0vZ2dnKyMhQXV2dCgsLW22LiooKVWkAAKCPC9kITVxcnE4//XRt2LBBkrR+/Xp5PB4VFBQoPT1dLtfR7OTxeOT3++X1ettsa0l9fb18Pl+zCQAAIGSBxuVy6c9//rMWLVqkiIgIPfnkk3r++edVUVGh+Pj4bw/odCouLk6VlZVttrUkLy9PsbGxwSk5OTlU5QMAAIuFLNAcPnxY11xzje677z55vV5dffXVuvbaayVJxphm6xpj5HA4Ttp2vCVLlqi2tjY4lZWVhap8AABgsZBdQ/Pf//3famho0C9+8QtJ0qOPPqqoqChNmzZN27dvD67n9/tVU1OjxMREJSUltdrWkoiICEVERISqZAAA0EeEbISmX79+8vv9wXljjAKBgC666CIVFRWpqalJklRcXCyXy6W0tDRNnTq11TYAAID2ClmgOf/88+Xz+bRs2TKVl5fr7rvvVnJyss4//3wlJCQoOztb5eXlys3N1axZsxQZGamMjIxW2wAAANorZIEmISFB69at03PPPafRo0dr8+bNeumllxQREaF169YpPz9fqampOnLkiJYuXXr04E5nq20AAADtFbJraCRp+vTpmj59+gnLJ06cqC1btrS4TVttAAAA7cG7nAAAgPUINAAAwHoEGgAAYL2QXkMDwF4LVxaFuwQA6DRGaAAAgPUYoTmVrJ7Tue2uXBvaOgAACDFGaAAAgPUINAAAwHoEGgAAYD0CDQAAsB6BBgAAWI9AAwAArEegAQAA1usbz6H58wIpsn+4qwAAAGHCCA0AALAegQYAAFiPQAMAAKxHoAEAANYj0AAAAOsRaAAAgPUINAAAwHoEGgAAYD0CDQAAsB6BBgAAWI9AAwAArEegAQAA1iPQAAAA6xFoAACA9Qg0AADAeq5wFwAAaNvClUUh3+ezCyaHfJ9AODFCAwAArEegAQAA1iPQAAAA6xFoAACA9Qg0AADAegQaAABgPW7bBizUHbfxAoDNGKEBAADWI9AAAADrhTTQNDY2atGiRYqOjtbYsWNVWFgoSfJ6vfJ4PHK73crKylJVVVVwm7baAAAA2iOkgeZXv/qVvvzySxUXF2vu3Lm6+uqrFQgENHv2bM2YMUOlpaVyu91avHixJLXZBgAA0F4OY4wJ1c5SU1P14osvasKECTp48KDWr1+vuLg4zZo1S/v375fL5ZLX61VGRoaqq6tVWFiomTNnttgWFRV10uP5fD7Fxsaq9vdXKCayf6g+Bo535dpwV4DjcFEwuop3OSGcgt/ftbWKiYkJyT5DNkJTUVGhL774Qps2bVJsbKwyMzM1YcIEvffee0pPT5fLdfSGKo/HI7/fL6/Xq4KCglbbWlJfXy+fz9dsAgAACFmg2bt3r5xOp95//31t2bJFY8aM0Q033KCKigrFx8d/e0CnU3FxcaqsrGyzrSV5eXmKjY0NTsnJyaEqHwAAWCxkgaaurk5+v1/Z2dkaNWqUbrnlFm3cuFGBQEDHn9UyxsjhcAR/bq3teEuWLFFtbW1wKisrC1X5AADAYiF7sF5sbKwkaciQIZKkoUOHyhij0047TTt27Aiu5/f7VVNTo8TERCUlJWn79u0ttrUkIiJCERERoSoZAAD0ESEboUlNTVX//v312WefSZIqKyvVr18/TZ06VUVFRWpqapIkFRcXy+VyKS0trc02AACA9gpZoHG73br88suVk5Ojzz//XL/+9a/1gx/8QBkZGUpISFB2drbKy8uVm5urWbNmKTIyss02AACA9grpc2iWL18uY4zOPfdcVVZW6sknn5TT6dS6deuUn5+v1NRUHTlyREuXLj168DbaAAAA2iukL6ccNmyYXn/99ROWT5w4UVu2bGlxm7baAAAA2oN3OQEAAOsRaAAAgPUINAAAwHoEGgAAYL2QXhSMPmr1nM5vy4stAQA9gBEaAABgPQINAACwHoEGAABYj0ADAACsR6ABAADWI9AAAADrEWgAAID1CDQAAMB6BBoAAGA9Ag0AALAegQYAAFiPQAMAAKxHoAEAANYj0AAAAOsRaAAAgPUINAAAwHoEGgAAYD0CDQAAsJ4r3AUAAHrewpVFId3fswsmh3R/QEcRaIAeEOovDwBAc5xyAgAA1iPQAAAA6xFoAACA9Qg0AADAegQaAABgPQINAACwHoEGAABYj0ADAACsR6ABAADWI9AAAADrEWgAAID1CDQAAMB6BBoAAGA9Ag0AALBeyAPN22+/LYfDobfeekuS5PV65fF45Ha7lZWVpaqqquC6bbUBAAC0V0gDTWNjo2666abgfCAQ0OzZszVjxgyVlpbK7XZr8eLFJ20DAADoCFcod7Zs2TINGzZM5eXlkqRNmzbpwIEDysnJkcvlUnZ2tjIyMlRXV6fCwsJW26KiokJZFgAA6ONCNkJTXl6uhx56SE899VRwWUFBgdLT0+VyHc1NHo9Hfr9fXq+3zTYAAICOCNkIza233qpFixZpzJgxwWUVFRWKj48PzjudTsXFxamysrLNttbU19ervr4+OO/z+UJVPgAAsFhIRmjWr1+vDz/8UPfcc88JbcaYE+YdDsdJ21qSl5en2NjY4JScnByC6gEAgO1CEmjWrl2rvXv3avjw4Ro8eLBqa2s1Y8YMJSUlqbq6Orie3+9XTU2NEhMT22xrzZIlS1RbWxucysrKQlE+AACwXEhOOT366KPKzs4Ozo8fP14rVqzQ8OHD9fDDD6upqUkul0vFxcVyuVxKS0uT3+9vta01ERERioiICEXJ6Cmr53R+2yvXhq4OAECfFpIRmvj4eI0aNSo4OZ1OJSYmKiMjQwkJCcrOzlZ5eblyc3M1a9YsRUZGttkGAADQEd36pGCn06l169YpPz9fqampOnLkiJYuXXrSNgAAgI4I6XNojqmpqQn+PHHiRG3ZsqXF9dpqAwAAaC/e5QQAAKxHoAEAANYj0AAAAOsRaAAAgPUINAAAwHoEGgAAYL1uuW0bAHBqWbiyKOT7fHbB5JDvE30XIzQAAMB6BBoAAGA9Ag0AALAegQYAAFiPQAMAAKxHoAEAANYj0AAAAOvxHBrgON3xPA0AQPdihAYAAFiPQAMAAKxHoAEAANYj0AAAAOsRaAAAgPUINAAAwHoEGgAAYD0CDQAAsB6BBgAAWI9AAwAArEegAQAA1iPQAAAA6/FySvReq+d0ftsr14auDgBAr8cIDQAAsB6BBgAAWI9AAwAArEegAQAA1iPQAAAA6xFoAACA9Qg0AADAegQaAABgPR6sB6uVlNW0uPw3K4t6thAAQFgxQgMAAKxHoAEAANYj0AAAAOuFNNB88cUXyszMVHR0tC688ELt3r1bkuT1euXxeOR2u5WVlaWqqqrgNm21AQAAtEdIA83111+vESNGaOvWrRo6dKhuvvlmBQIBzZ49WzNmzFBpaancbrcWL14sSW22AQAAtFfI7nJqaGjQm2++qa1bt2rkyJFauHCh5s2bp02bNunAgQPKycmRy+VSdna2MjIyVFdXp8LCwlbboqKiQlUaAADo40I2QtPY2KhHHnlEKSkpkqT9+/fL7XaroKBA6enpcrmOZiePxyO/3y+v19tmW0vq6+vl8/maTQAAACELNFFRUfr5z38ut9utxsZGPfHEE5o/f74qKioUHx//7QGdTsXFxamysrLNtpbk5eUpNjY2OCUnJ4eqfAAAYLGQ3+XU1NSkq666Sk6nU7m5uZIkY0yzdYwxcjgcJ2073pIlS1RbWxucysrKQl0+AACwUEifFBwIBDR37lx98cUX+p//+R+53W4lJSVp+/btwXX8fr9qamqUmJjYZltLIiIiFBEREcqSAQBAHxDSEZrc3Fzt3LlTb775poYMGSJJmjp1qoqKitTU1CRJKi4ulsvlUlpaWpttAAAA7RWyQFNRUaFly5bpt7/9rSSppqZGNTU1ysjIUEJCgrKzs1VeXq7c3FzNmjVLkZGRbbYBAAC0V8gCzYYNG+Tz+TRlyhTFxcUFpz179mjdunXKz89Xamqqjhw5oqVLlx49uNPZahsAAEB7OczxV+VaxOfzKTY2VrW/v0Ixkf3DXQ7CoNW3bX/ngZ4tBEDIPbtgcrhLQDcJfn/X1iomJiYk++RdTgAAwHoEGgAAYD0CDQAAsF5In0MD9Bb/Vvmfnd6W628AwD6M0AAAAOsRaAAAgPU45QQA6JUWriwK+T65FbzvYoQGAABYj0ADAACsR6ABAADWI9AAAADrEWgAAID1uMsJPaq1l0kCANAVjNAAAADrEWgAAID1CDQAAMB6XEMDHIcXWwKAfRihAQAA1iPQAAAA6xFoAACA9Qg0AADAegQaAABgPQINAACwHoEGAABYj0ADAACsR6ABAADWI9AAAADr8eoDIIR4bQIAhAeBBq0qKasJdwkAEFILVxaFdH/PLpgc0v2h8zjlBAAArEegAQAA1iPQAAAA6xFoAACA9bgouA/hIl67dfYOKe6OAgBGaAAAQB9AoAEAANYj0AAAAOsRaAAAgPUINAAAwHphv8vJ6/Xq2muv1Y4dOzR16lS98MILGjZsWLjL6nbckQQAQOg4jDEmXAcPBAI688wzddVVV+mGG27QzTffrEGDBmnVqlXt2t7n8yk2Nla1v79CMZH9u7na0CLQwHbcLg50j1Ph/VDB7+/aWsXExIRkn2Edodm0aZMOHDignJwcuVwuZWdnKyMjQ3V1dYqKigpnaQAAwCJhDTQFBQVKT0+Xy3W0DI/HI7/fL6/Xq6lTp56wfn19verr64PztbW1kiTf4cZurfOj/6vt1v0DNvrp7jvDctynh3XuAYSALXw+X7hL6HbHPmMoTxKFNdBUVFQoPj4+OO90OhUXF6fKysoW18/Ly9N99913wvLkf8/vthoB9DZvhrsAoFu9cFO4K+g5+/fvV2xsbEj2FfaLgo9PZ8YYORyOFtddsmSJbrvttuB8TU2NRo4cqT179oTsF4LO8fl8Sk5OVllZWcjOh6Jz6Iveg77oXeiP3qO2tlYjRozQkCFDQrbPsAaapKQkbd++PTjv9/tVU1OjxMTEFtePiIhQRETECctjY2P5j7OXiImJoS96Cfqi96Avehf6o/dwOkP39JiwPodm6tSpKioqUlNTkySpuLhYLpdLaWlp4SwLAABYJqyBJiMjQwkJCcrOzlZ5eblyc3M1a9YsRUZGhrMsAABgmbAGGqfTqXXr1ik/P1+pqak6cuSIli5d2u7tIyIilJ2d3eJpKPQs+qL3oC96D/qid6E/eo/u6IuwPlgPAAAgFHiXEwAAsB6BBgAAWI9AAwAArEegAQAA1uvVgcbr9crj8cjtdisrK0tVVVUn3WbDhg06++yzFRUVpblz5+rQoUM9UGnf15m++PDDDzVx4kTFxMToRz/6kQ4cONADlfZ9nemLY55//nk5HA59+eWX3VfgKaQzfeHz+TR79mxFRUXp/PPP186dO3ug0r6vM32xadMmjR8/XjExMfrJT36ir7/+ugcqPXXs3btXmZmZKikpadf6Xf3+7rWBJhAIaPbs2ZoxY4ZKS0vldru1ePHiNrepqanRT37yE912223atm2bvvzySz344IM9VHHf1Zm+CAQCuvLKK5WVlaWPPvpIVVVVuvfee3uo4r6rM31xTE1Nje64445urvDU0dm+uOOOOxQdHa1t27ZpwoQJuummU+jFPd2kM33xzTffaObMmbr99tu1bds2HTlyhL+PEFq0aJGGDx+ut99+u13rh+T72/RSb775pomJiTGNjY3GGGM+/PBD43a7zcGDB1vd5rnnnjNjx44Nzr/44otmxIgR3V5rX9eZvigtLTWSgus89dRTZvz48T1Sb1/Wmb445sYbbzQzZ840ksyuXbu6udK+rzN9cejQIRMXF2f27dtnjDGmoqLCvPrqqz1Sb1/Wmb4oLCw0ERERwfnVq1c3+/5A11RXV5tdu3YZSaa4uPik64fi+7vXjtAUFBQoPT1dLtfR1015PB75/X55vd42t5kyZUpw/vzzz9eePXtUVlbW7fX2ZZ3pi8jISD3++OOKioqSdPSNqm63u0fq7cs60xfS0dN/69at0yOPPNITZZ4SOtMXW7ZsUXR0tB577DENGjRI8+bNU3p6ek+V3Gd1pi/OOOMMuVwubdq0SX6/X6+//ro8Hk8PVdz3xcfHa9SoUe1ePxTf37020FRUVCg+Pj4473Q6FRcXp8rKynZvM3ToUElqcxucXGf6Yvjw4br11lslHX2r6jPPPKP58+d3e619XWf6IhAI6MYbb9QDDzyghISEnijzlNCZvti7d6+qq6tVV1enjz/+WA6HQ0uWLOmJcvu0zvTF0KFDtXLlSl100UWKjIzUW2+9pSeeeKInykULQvH93WsDjSSZ4x5ibIyRw+Fo9zbHfj7ZNji5zvSFJNXV1emyyy5TWlqabrzxxu4q75TS0b54+umn5XA4dN1113V3aaecjvZFXV2dGhoa9NBDDyklJUXXXXed1q9f391lnhI62hd79+7VLbfcoueee04ffPCBJk+e3O7r0dA9uvr93WsDTVJSkqqrq4Pzfr9fNTU1SkxMbPc2+/fvl6Q2t8HJdaYvJOnQoUOaPn26IiMj9V//9V8hfU38qaozffHnP/9ZH330kYYMGaKRI0dKksaPH6/Nmzd3e719WWf6IjY2VgMHDgyefh06dGjw3yl0Xmf/LkaPHq0FCxbo3HPP1UMPPaQXXnhBtbW1PVEyjhOK7+9e+w0zdepUFRUVqampSZJUXFwsl8ultLS0Nrd59913g/P/+7//q1GjRum0007r9nr7ss70hSTddNNNioqK0iuvvKKBAwf2RKl9Xmf6Ys2aNdq+fbtKSkr0zjvvSJJee+01TZo0qUdq7qs60xdjxoxRXV2dvvrqK0lHh9P5H66u60xf9OvXT36/Pzjf1NQkY0zwOhz0rJB8f3foEuIe5Pf7zRlnnGHuuusuU1ZWZi677DJz9dVXG2OMqa2tNQ0NDSdsU1NTY2JjY83y5cvNrl27THp6uvnP//zPni69z+lMX5SUlJjY2Fjz2Wefma+//jo4oWs60xf/6Ouvv+YupxDpbF+cd9555pprrjGff/65ueCCC8wtt9zSk2X3SZ3pi61bt5r+/fubNWvWmN27d5t58+aZKVOm9HTpfZ6Ou8upO7+/e22gMeborXfjx483ERER5pJLLjFVVVXGGGNGjhxpXn755Ra32bBhg0lNTTVut9vMmTPH1NXV9WDFfVdH+yInJ8dIOmFC13Xm7+IYAk1odaYvduzYYSZPnmyioqLMjBkzzIEDB3qw4r6rM33xpz/9yaSmppqoqChz8cUXm88//7wHKz41HB9ouvP72/H/DwgAAGCtXnsNDQAAQHsRaAAAgPUINAAAwHoEGgAAYD0CDQAAsB6BBgAAWI9AAwAArEegAQAA1iPQAAAA6xFoAACA9f4fr96wWEAhpmMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "class GibbsSampling:\n",
    "    def __init__(self, target_dist, j, m=1e4, n=1e5):\n",
    "        \"\"\"\n",
    "        Gibbs Sampling 算法\n",
    "\n",
    "        :param target_dist: 目标分布\n",
    "        :param j: 变量维度\n",
    "        :param m: 收敛步数\n",
    "        :param n: 迭代步数\n",
    "        \"\"\"\n",
    "        self.target_dist = target_dist\n",
    "        self.j = j\n",
    "        self.m = int(m)\n",
    "        self.n = int(n)\n",
    "\n",
    "    def solve(self):\n",
    "        \"\"\"\n",
    "        Gibbs Sampling 算法求解\n",
    "        \"\"\"\n",
    "        # (1) 初始化\n",
    "        all_samples = np.zeros((self.n, self.j))\n",
    "        # 任意选择一个初始值\n",
    "        x_0 = np.random.random(self.j)\n",
    "        # x_0 = np.array([0.5, 0.5])\n",
    "        # (2) 循环执行\n",
    "        for i in range(self.n):\n",
    "            x = x_0 if i == 0 else all_samples[i - 1]\n",
    "            # 满条件分布抽取\n",
    "            for k in range(self.j):\n",
    "                x[k] = self.target_dist.sample(x, k)\n",
    "            all_samples[i] = x\n",
    "        # (3) 得到样本集合\n",
    "        samples = all_samples[self.m:]\n",
    "        # (4) 计算函数样本均值\n",
    "        dist_mean = samples.mean(0)\n",
    "        dist_var = samples.var(0)\n",
    "        return samples[self.m:], dist_mean, dist_var\n",
    "\n",
    "    @staticmethod\n",
    "    def visualize(samples, bins=50):\n",
    "        \"\"\"\n",
    "        可视化展示\n",
    "        :param samples: 抽取的随机样本集合\n",
    "        :param bins: 频率直方图的分组个数\n",
    "        \"\"\"\n",
    "        fig, ax = plt.subplots()\n",
    "        ax.set_title('Gibbs Sampling')\n",
    "        ax.hist(samples[:, 0], bins, alpha=0.7, label='$\\\\theta$')\n",
    "        ax.hist(samples[:, 1], bins, alpha=0.7, label='$\\\\eta$')\n",
    "        ax.set_xlim(0, 1)\n",
    "        ax.legend()\n",
    "        plt.show()\n",
    "\n",
    "        \n",
    "        \n",
    "class TargetDistribution:\n",
    "    \"\"\"\n",
    "    目标概率分布\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self):\n",
    "        # 联合概率值过小，可对建议分布进行放缩\n",
    "        self.c = self.__select_prob_scaler()\n",
    "\n",
    "    def sample(self, x, k=0):\n",
    "        \"\"\"\n",
    "        使用接受-拒绝方法从满条件分布中抽取新的分量 x_k\n",
    "        \"\"\"\n",
    "        theta, eta = x\n",
    "        if k == 0:\n",
    "            while True:\n",
    "                new_theta = np.random.uniform(0, 1 - eta)\n",
    "                alpha = np.random.uniform()\n",
    "                if (alpha * self.c) < self.__prob([new_theta, eta]):\n",
    "                    return new_theta\n",
    "        elif k == 1:\n",
    "            while True:\n",
    "                new_eta = np.random.uniform(0, 1 - theta)\n",
    "                alpha = np.random.uniform()\n",
    "                if (alpha * self.c) < self.__prob([theta, new_eta]):\n",
    "                    return new_eta\n",
    "\n",
    "    def __select_prob_scaler(self):\n",
    "        \"\"\"\n",
    "        选择合适的建议分布放缩尺度\n",
    "        \"\"\"\n",
    "        prob_list = []\n",
    "        step = 1e-3\n",
    "        for theta in np.arange(step, 1, step):\n",
    "            for eta in np.arange(step, 1 - theta + step, step):\n",
    "                prob = self.__prob((theta, eta))\n",
    "                prob_list.append(prob)\n",
    "        searched_max_prob = max(prob_list)\n",
    "        upper_bound_prob = searched_max_prob * 10\n",
    "        return upper_bound_prob\n",
    "\n",
    "    @staticmethod\n",
    "    def __prob(x):\n",
    "        \"\"\"\n",
    "        P(X = x) 的概率\n",
    "        \"\"\"\n",
    "        theta = x[0]\n",
    "        eta = x[1]\n",
    "        p1 = (theta / 4 + 1 / 8) ** 14\n",
    "        p2 = theta / 4\n",
    "        p3 = eta / 4\n",
    "        p4 = (eta / 4 + 3 / 8)\n",
    "        p5 = 1 / 2 * (1 - theta - eta) ** 5\n",
    "        p = (p1 * p2 * p3 * p4 * p5)\n",
    "        return p\n",
    "\n",
    "# 收敛步数\n",
    "m = 1e3\n",
    "# 迭代步数\n",
    "n = 1e4\n",
    "\n",
    "# 目标分布\n",
    "target_dist = TargetDistribution()\n",
    "\n",
    "# 使用 Gibbs Sampling 算法进行求解\n",
    "gibbs_sampling = GibbsSampling(target_dist, 2, m, n)\n",
    "\n",
    "samples, dist_mean, dist_var = gibbs_sampling.solve()\n",
    "\n",
    "print(f'theta均值：{dist_mean[0]}, theta方差：{dist_var[0]}')\n",
    "print(f'eta均值：{dist_mean[1]}, eta方差：{dist_var[1]}')\n",
    "\n",
    "# 对结果进行可视化\n",
    "GibbsSampling.visualize(samples, bins=20)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a3c20f5",
   "metadata": {},
   "source": [
    "# 平稳分布案例"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "664fd532",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.23076935 0.30769244 0.46153864]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAGbCAYAAAABeQD9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABGFElEQVR4nO3deXxU9d3+/9csyWQme0gghASSEKxiwYCIIkEpKvUruFSsilv1VlqBu7dL1f7wrqLcKrW2VlrcaluXSrXSilpwqVWwAi5o4oKgsiQkgZCEkD1kmzm/PyYZAmSZhEnOTLiej57HzFnnPSfHzsU5n/M5FsMwDERERESCjNXsAkREREQ6o5AiIiIiQUkhRURERIKSQoqIiIgEJYUUERERCUoKKSIiIhKUFFJEREQkKCmkiIiISFCym13A0fB4POzZs4fo6GgsFovZ5YiIiIgfDMOgtraWlJQUrNauz5eEdEjZs2cPaWlpZpchIiIifVBUVERqamqX80M6pERHRwPeLxkTE2NyNSIiIuKPmpoa0tLSfL/jXQnpkNJ+iScmJkYhRUREJMT01FRDDWdFREQkKCmkiIiISFBSSBEREZGgpJAiIiIiQUkhRURERIKSQoqIiIgEJYUUERERCUoKKSIiIhKUFFJEREQkKCmkiIiISFBSSBEREZGgpJAiIiIiQSmkHzAoEmwMw8BjeGg1Wmn1HDa0TXN73LR4Wg4ZdxtuPIany8FtuPHgwePxeF87Wcb3+XgwDAMDo8tXj+HxjQOHzGsfb3tL25q+7Xdcz/e9OTjecd4h2zs48chph+3DI6Z1smxX6/ujs88Qkc4tzF5IVHiUKZ+tkCLHFI/hoba5lprmGu/Q5H2tba6lsbWRRncjTe4mmlqbaHQ3HpzW2kSTu8n3/vDlWjwt3uDhaTX7K4qIBNT1464nCoUUkV6raa5hT90edtftpuJAxRHh4/D3dc11R/Uv8L6yW+zYrd7BZrUdMm632rFarNgsNiwWi/cV76vVYu10OHxZi8WCBQtWi9X36HOrxeqd1z6/wzKHzG9b3kLba4fxw6f5pluOnNc+frjO5ne23UPW6WxbnU7q/jHv3c3v6RHx/a2n2kWChdPuNO2zFVIkqNU217Knbg/FdcXsqdvjCyTt72tbavu0XafdSXR4NDHhMcQ6YokOj8Zpc+KwO3DYHETYInDY215tDiLs3tfOpkXYIwi3hh8MHYcFkPZAISIivaOQIqZrcbfwVcVXfFXxFcW1bWGk3htGapt7DiEJEQmMiBpBojORmPAYYhwx3tcu3seGxxJmCxuAbyYiIkdDIUUGXF1zHZ+Vf0ZuaS65Zbls3reZJndTl8snRCSQEplCSlQKI6JG+F5HRI0gOTIZV5hrAKsXEZGBopAi/a68oZxPyz4lrzSP3LJcvq381nc3SruEiATGJ40nIyaDlKiDgWR45HCFEBGRY5RCigSUYRgU1BT4zpLkluZSXFd8xHKpUalMHDaRiUMnMnHYRNJj0tVuQ0REDqGQIgHxWdln/PXrv/JRyUfsb9x/yDwLFr6T8B0mDp3IhGETmDh0IkNdQ02qVEREQoVCivRZq6eVdwrf4bktz/FF+Re+6eHWcMYljfOdJTkp6SSiw6NNrFREREKRQor0Wl1zHS9ve5kVW1ewp34PAGHWMGZnzuYHY37AiUNOJNwWbnKVIiIS6hRSxG8ldSWs2LqCf2z7B3UtdQDEO+K57PjLuOw7l5HoTDS5QhERGUwUUqRHX5Z/yXNbnuPtXW/jNtwAZMRmcM3Ya5idOZsIe4TJFYqIyGCkkCKdcnvcrCtax3NbniO3LNc3/dThp3LN2GvIGZHj615dRESkP/TqVyY3N5fs7GycTiczZ86krKzM73Wfe+45LBYLBQUFvmnJyckHnylisTBp0qTelCP9oKGlgRVbVzB71WxuXnczuWW52K12Lhh9AX8//+/8ceYfOSP1DAUUERHpd36fSfF4PMyZM4crr7yS1atXs3DhQm655RZWrFjR47pVVVXcfvvtnU7funUrycnJ3mLsOrFjppe+eYlHch/xdUUf64jl0uMuZe7xc0lyJZlcnYiIHGsshmH49UjYtWvXctFFF1FRUYHdbic3N5ecnBzKy8uJjIzsdt0FCxZQUlLCK6+8Qn5+Punp6TQ1NREREUFTUxPh4X27E6SmpobY2Fiqq6uJiYnp0zbEeyvxLz/+JX/75m8AjIoZxdUnXM35o89Xb68iIhJw/v5++33OfsOGDUyePNl3tiM7Oxu3201ubm6363366aesXLmSX/3qV4dMr6qqIiIigrlz5+J0OjnjjDPYvXt3t9tqamqipqbmkEGOTk1zDQvfWcjfvvkbFizcPPFmXrvoNS47/jIFFBERMZXfIWXv3r0kJh68xdRqtRIfH09paWmX63g8HubPn899991HUtKhlwuqqqpobGzkrLPOYsuWLVitVm677bZua1i6dCmxsbG+IS0tzd/ypRNFNUVc9fpVbNyzEafdyW+/91uuH3e92puIiEhQ6NWv0eFXhgzD6PZ5K08++SQWi4V58+YdMS8zM5M9e/awYMECMjIymD9/PmvXru328xctWkR1dbVvKCoq6k350sGnpZ9yxetXkF+dz1DXUJ4991nOGnmW2WWJiIj4+B1Shg8fTnl5uW/c7XZTVVXla/TamZdeeokvvviChIQERo0aBcD48eNZv349YWFhDB8+3LdsfHx8j5dvHA4HMTExhwzSe69sf4Ub/nUDVU1VnDjkRF6Y9QInDDnB7LJEREQO4XdImTZtGps2baK1tRWAvLw87HY7EyZM6HKdF154ga1bt/LZZ5/x/vvvA/D6668zadIkHn30UWbMmOFbtrCwkPT09D5+DfGHx/DwyKePcNeGu2j1tHLOqHN4+tyn9bA/EREJSn6HlJycHJKSkli8eDHFxcUsWbKEiy++GJfLRU1NDS0tLUesk5ycTHp6Ounp6YwcORKA1NRUIiIimD59Ohs3bmT16tV88803LFu2jGuvvTZgX0wO1dDSwK3rbuVPm/8EwI/H/5hfn/lrnHanyZWJiIh0zu+QYrVaWblyJatXryYrK4vGxkYefvhhwHsJZ82aNb364BNPPJEnnniChQsXMm3aNM4++2xuueWW3lUvfimtL+XaN6/lncJ3CLOG8UDOA/x0wk/VQFZERIKa3/2kBCP1k9Kzr/Z9xU/f/SnlB8pJiEjgke89woShXV+iExER6W/+/n6ri9dB7F8F/+J/1/8vje5GsuKyWH7WckZEjTC7LBEREb8opAxChmHw1JdP8fu83wOQMyKHh854iKjwKJMrExER8Z9CyiDT7G5m8cbFrN65GoCrTriKn036GXar/tQiIhJa9Ms1iDS0NHDjv28krywPm8XGnafeyaXfudTsskRERPpEIWUQWf7ZcvLK8ogOi+Y303/DlJQpZpckIiLSZ7oHdZDYUrGFFVtXAPDgGQ8qoIiISMhTSBkEWj2t3PvBvXgMD+emn8u01GlmlyQiInLUFFIGgRe+foEtFVuIDovm55N/bnY5IiIiAaGQEuL21u/13Wp888k3k+hMNLkiERGRwFBICWGGYXD/R/dzoPUAE4ZO4JLjLjG7JBERkYBRSAlh7xS+w7qidditdu4+7W49i0dERAYV/aqFqLrmOpZ+tBSA6068jqz4LJMrEhERCSyFlBD1u7zfUXagjJHRI/nx+B+bXY6IiEjAKaSEoC/Kv+DFr18E4K4pdxFhjzC5IhERkcBTSAkxLZ4W7v3gXgwMzs88n9OGn2Z2SSIiIv1CISXEPL/leb6t/JZYRyy3nXKb2eWIiIj0G4WUEFJcW8xjnz0GwM9O/hkJEQkmVyQiItJ/FFJChGEY3PfRfTS6Gzkl+RQuyrrI7JJERET6lUJKiHir4C027N5AmDWMu067C4vFYnZJIiIi/UohJQRUN1Xzy49/CcC8cfPIiM0wuSIREZH+p5ASApblLqOisYKM2AyuH3e92eWIiIgMCIWUIJdXlsfKb1cCcPdpdxNuCze5IhERkYGhkBLEWtwt3LvxXgB+kPUDJiVPMrkiERGRgaOQEsSe/uppdlTvICEigZ9N+pnZ5YiIiAwohZQgVVhTyJOfPwnA7afcTqwj1uSKREREBpZCShAyDIMlHy6h2dPMlOFTmJUxy+ySREREBpxCShBavXM1H5V8hMPmUJ8oIiJyzFJICTJVjVU8tOkhAG486UbSYtJMrkhERMQcCilBZlneMiqbKsmKy+JHJ/7I7HJERERMo5ASRA60HmDNzjUA3HnqnYRZw0yuSERExDwKKUFk/e71HGg9QEpkCpOGqU8UERE5timkBJG3C94G4JxR56ixrIiIHPMUUoJEY2sj7xW/B8DM9JkmVyMiImI+hZQgsWHPBhpaG0iOTGZc4jizyxERETGdQkqQeHuXLvWIiIh0pJASBJrcTawrWgfAzFG61CMiIgIKKUHhgz0fUN9Sz1DXUMYnjTe7HBERkaCgkBIEOl7qsVr0JxEREQGFFNM1u5tZW7gW8IYUERER8VJIMdmHJR9S21JLojOR7KRss8sREREJGgopJmu/1HP2yLOxWW0mVyMiIhI8FFJM1OJp4d3CdwF14CYiInI4hRQTfVzyMTXNNSREJDBx6ESzyxEREQkqCikm0qUeERGRrimkmKTF08I7he8AcE667uoRERE5nEKKST7Z+wlVTVXEO+KZNGyS2eWIiIgEHYUUk/xr178AmDFyBnar3eRqREREgo9CiglaPa0H7+rRs3pEREQ6pZBigtzSXPY37ifWEcspw08xuxwREZGgpJBiAt+lnrQZhFnDTK5GREQkOCmkDDC3x82/d/0bUAduIiIi3VFIGWB5ZXlUNFYQHR7Nqcmnml2OiIhI0FJIGWDtl3q+l/Y9wmy61CMiItIVhZQB5DE8vks930//vsnViIiIBDeFlAH0efnnlB8oJyositOGn2Z2OSIiIkFNIWUA/avAe6lnetp0wm3hJlcjIiIS3BRSBojH8PgeKKgO3ERERHqmkDJAvij/gtKGUiLDIjl9xOlmlyMiIhL0FFIGSPtZlDNTz8Rhc5hcjYiISPBTSBkAhmHoUo+IiEgvKaQMgM37NlNSX4LT7mTqiKlmlyMiIhISehVScnNzyc7Oxul0MnPmTMrKyvxe97nnnsNisVBQUOCb9tZbb3HccccRGRnJ5ZdfTkNDQ2/KCRkdL/VE2CNMrkZERCQ0+B1SPB4Pc+bMYfbs2Wzbtg2n08ktt9zi17pVVVXcfvvtR0y79NJLufXWW/nqq68oKCjg/vvv7131IcAwDF8vs+eMOsfkakREREKH3yHlvffeY//+/dxzzz2kpqayePFiVq1aRX19fY/r3nnnnZx++qF3tKxatYrU1FRuvPFG0tPTueOOO3j++ed7/w2C3Jb9W9hdt5sIWwQ5I3LMLkdERCRk+B1SNmzYwOTJk7Hb7QBkZ2fjdrvJzc3tdr1PP/2UlStX8qtf/eqI7XUMLqeeeiqFhYUUFRX1pv6g93aB91LPtNRpuMJcJlcjIiISOvwOKXv37iUxMfHgilYr8fHxlJaWdrmOx+Nh/vz53HfffSQlJXW7vSFDhgB0u72mpiZqamoOGYKZ7uoRERHpu141nDUM44hxi8XS5fJPPvkkFouFefPm9bi99vfdbW/p0qXExsb6hrS0tN6UP+C+qfyGwtpCHDYHZ6SeYXY5IiIiIcXvkDJ8+HDKy8t94263m6qqKpKTk7tc56WXXuKLL74gISGBUaNGATB+/HjWr19/xPYqKioAut3eokWLqK6u9g3Bfmmo/Vk9OSNydKlHRESkl+z+Ljht2jQefPBBWltbsdvt5OXlYbfbmTBhQpfrvPDCCzQ2NgJQU1PDSSedxOuvv86kSZPYuXMnS5cu9S374Ycfkp6ezogRI7rcnsPhwOEIjd5aO17q0V09IiIivef3mZScnBySkpJYvHgxxcXFLFmyhIsvvhiXy0VNTQ0tLS1HrJOcnEx6ejrp6emMHDkSgNTUVCIiIrjwwgspKSnh8ccfp6CggIceeoirrroqcN/MZNuqtlFQU0C4NZwzU880uxwREZGQ43dIsVqtrFy5ktWrV5OVlUVjYyMPP/ww4L2Es2bNml59cGxsLC+99BIPP/wwY8eOJSMjg0WLFvWu+iDWfhbl9BGnExUeZXI1IiIiocfvyz0AEydO5PPPPz9iesdeZLsSFxd3RMPbmTNnsm3btt6UEDLa26Porh4REZG+0bN7+sGOqh3srN6J3Wpnetp0s8sREREJSQop/aC9G/zTU04nOjza5GpERERCk0JKP9ClHhERkaOnkBJg+dX5bK/ajt2iSz0iIiJHQyElwDbt3QTAycknE+uINbkaERGR0KWQEmDbKr13K41NGGtyJSIiIqFNISXAdlTvAGB03GiTKxEREQltCimH8XgM3txcwqNrt9PY4u71+juqvCElKy4r0KWJiIgcU3rVmduxwGKBn//jS6oPtHD2CcP4TrL/txBXHKhgf+N+LFjIiM3oxypFREQGP51JOYzFYiE9MRKA/H11vVq3/SzKiKgReuqxiIjIUVJI6URmW0jZua++V+ttq/I2ms2K16UeERGRo6WQ0omM9jMp5b0LKWqPIiIiEjgKKZ3whZRenklpDym6s0dEROToKaR0oi8hxTAM3+WeMXFj+qUuERGRY4lCSifaQ0pFfTPVDS1+rVN+oJza5lqsFivpsen9WJ2IiMixQSGlE5EOO8NiHADkV/h3NmV75XYARkaPxGFz9FttIiIixwqFlC5k9PI25O1V3pCiRrMiIiKBoZDShYzEKMD/O3zUHb6IiEhgKaR0obd9pbRf7lEfKSIiIoGhkNKF3tzhYxjGwcs9sQopIiIigaCQ0oWMpIMhxTCMbpctqS+hobUBu9XOqJhRA1GeiIjIoKeQ0oW0eBc2q4WGZjdltU3dLtt+FiU9Jp0wW9hAlCciIjLoKaR0IdxuJS3eCcDOHhrPqqdZERGRwFNI6Ya/7VJ0+7GIiEjgKaR0w3cbcg99pSikiIiIBJ5CSjc6Np7tisfwsLNqJ6CQIiIiEkgKKd3wp6+U3bW7aXQ3Em4NJy06baBKExERGfQUUrrR3ialsKKBVren02XaL/VkxGZgs9oGrDYREZHBTiGlG8kxEUSEWWn1GBRXHuh0GV97FPU0KyIiElAKKd2wWi2kD+m+XYoazYqIiPQPhZQeZCZ13y5FIUVERKR/KKT04GBfKUfehtzqaSW/Oh9QR24iIiKBppDSg4N9pRx5JqWotogWTwtOu5MRUSMGujQREZFBTSGlB74zKZ10jd9+qSczNhOrRbtSREQkkPTL2oP2vlL2VDdyoNl9yLz2kKJLPSIiIoGnkNKD+Mhw4lzeJxsXVBx6NqX9wYJj4sYMeF0iIiKDnUKKH7p60OD2Sp1JERER6S8KKX7oLKS0uFvYVbML0O3HIiIi/UEhxQ8ZnXToVlBTQKvRSmRYJMmRyWaVJiIiMmgppPihs6cht7dHGR03GovFYkpdIiIig5lCih86u9zTfmePGs2KiIj0D7vZBYSC9uf37K9vpqqhmThXuG4/FhExmcfjobm52ewypBNhYWHYbLaj3o5Cih8iHXaSYyLYW9NI/r56JowMP+Ryj4iIDKzm5mby8/PxeDxmlyJdiIuLIzk5+aiaRCik+CkjMdIXUsaOcFFYWwjoco+IyEAzDIOSkhJsNhtpaWlYrWq5EEwMw6ChoYGysjIAhg8f3udtKaT4KSMpkg92VpC/r5786jo8hoeY8BgSnYlmlyYickxpbW2loaGBlJQUXC6X2eVIJ5xOJwBlZWUMHTq0z5d+FFL81N49/s599Wyv2gN4+0fRnT0iIgPL7fY+oiQ8PNzkSqQ77QGypaWlzyFF58j81PFBg+09zaoTNxER8+gficEtEH8fhRQ/dbwNWXf2iIiI9D+FFD+lJbiwWS0caHHzzf5tAIyJV6NZERE5OuvWrSM9PX3AP7e2tpbZs2fjcrlITk5m06ZNfq9rsVgoKCjov+LaqE2Kn8JsVkYmuMivqGRvg7dNis6kiIhIqHrmmWcoKSlh+/btVFdXk5CQYHZJR1BI6YWMxEh21X0DQEJEAgkRwfcHFRER8UdFRQXjxo0jJSWFlJQUs8vplC739EJGYiRWx15AjWZFRIKFYRg0NLeaMhiG0ata8/LymDJlCtHR0UybNo2tW7f65r322muMGjWKhIQEli9f7pu+fft2zj77bGJiYjj55JMPuSzT3bx77rmHa6+9ljfffJPs7GzuuOMOAF588UUsFgv33nsvzz77LBaLheOPPx6AgoKCIxq8pqens27dul59z0DRmZReyEiMxObwdk6jSz0iIsHhQIubsXe/Zcpnb1nyfVzh/v2U1tTUcO6553LzzTfz97//nV//+tfMnTuXRx55hIqKCh588EFef/111q5dy89+9jNuuOEG7HY7F1xwAT//+c/585//zMsvv8zcuXPZtm0bbre7y3ntQWPz5s3ccccdLF68mOzsbADmzJlDZWUlv/zlLyksLOSxxx4LSBf2/UEhpRcydSZFRET6aPXq1SQkJLBo0SIAFi9eTE5ODgB1dXU88cQTnHjiiYwZM4af/vSnlJaWUlxczNdff81NN93k2051dTUlJSXk5+d3Oa/98s3mzZv5+uuvD2mYGxYWRlxcHBEREYSHhxMXF9f/X76PFFJ6ISMpEqujFID0mEyTqxEREQBnmI0tS75v2mf7q6io6JCwEBcXx5w5c1i3bh3x8fGMGzcOONhJnWEYFBcXM2rUKNauXXvItpKSknj//fe7nNfu/PPPP+o7h+rr649q/aOhNim9EOloxRpWDYCTESZXIyIi4L0d1hVuN2XoTYdlaWlph9y2W1dXx3e/+1327t1LTExMp+ukpqb6upZPT08nNTWVX/3qV1RUVHQ7r11kZGSv9yUc7NV3165d7Nu3r1fbCCSFlF7YWeN98rGnJYbyau06ERHx36xZs6ioqOCBBx6guLiY+++/H7fbzbBhw7pcZ/LkyYwYMYLbbruNoqIili5dypo1a0hKSup2Xl8NHTqUsLAw1q9fT2trK7feequpjx/QL20v7KhqCylNw9hZbt7pLxERCT2xsbG88cYbvPbaa5xwwgl88MEHrFq1qtuzMWFhYfzzn//k66+/ZuzYsaxevZpXX30Vm83W7by+cjqdPPjgg1x++eVMmDCB884776ieYny0LEZv758KIjU1NcTGxlJdXd3lqbJAevDjB3l+6/M0V+RwScYC7v/BuH7/TBEROVRjYyP5+flkZGQQERFhdjnShe7+Tv7+futMSi+0P7PH0zSM/H06kyIiItKfFFJ6of1yj1shRUREpN8ppPipuqma8gPlAHiah1FS3UhDc6vJVYmIiAxevQopubm5ZGdn43Q6mTlzJmVlZT2us2HDBsaNG0d0dDSzZs06ZJ3k5GQsFotvmDRpUu+/wQBpv9QzPHI48RFRABTsazCzJBERkUHN75Di8XiYM2cOs2fPZtu2bTidTm655ZZu12lsbOSHP/wht99+u+/5BO097QFUVVWxdetWKisrqaysNO3ZAP7YXukNKVlxWWQkeu871yUfERGR/uN3j7Pvvfce+/fv55577sFut/u6862vr++ys5j8/HwmTZrENddcA3jvEf/LX/4CQFNTE01NTWRmZpp6D7a/2s+kZMVl4UyMIrewivx9dSZXJSIiMnj5fSZlw4YNTJ48Gbvdm2uys7Nxu93k5uZ2uc4JJ5zAa6+9BkBtbS2vvPIKEyZMALxnUSIiIpg7dy5Op5MzzjiD3bt3H8136Vc7qr2NZkfHjSYzyRvKdupMioiISL/xO6Ts3buXxMTEgytarcTHx1NaWtrjuoWFhcTFxVFQUMD9998PeENKY2MjZ511Flu2bMFqtXLbbbd1u52mpiZqamoOGQaK73JPvC73iIiIDIReNZw9vN83wzD8em5BSkoK69evJzY2lrvuuguAzMxM9uzZw4IFC8jIyGD+/PlHPCTpcEuXLiU2NtY3pKWl9ab8Pqs4UEFlUyUWLGTGZpI+RCFFRESkv/kdUoYPH055eblv3O12U1VVRXJycpfrFBUVsXXrVux2O1OmTOHOO+9kxYoVgLer345d7cbHx/d4ZmTRokVUV1f7hqKiIn/LPyrt7VFSo1Nx2p2kJ7oAqGpoobK+eUBqEBEROdb4HVKmTZvGpk2baG319g2Sl5eH3W73tTHpzBtvvMH111/vG7dYLL5nCjz66KPMmDHDN6+wsLDHx0k7HA5iYmIOGQZCe0gZHTcaAFe4neGx3i5+1S5FRESkf/gdUnJyckhKSmLx4sUUFxezZMkSLr74YlwuFzU1NbS0tByxzowZM8jLy2PVqlXs2rWLRx99lPPOOw+A6dOns3HjRlavXs0333zDsmXLuPbaawP2xQKpvafZrLgs3zS1SxERCRKGAc315gy9fPxdXl4ep5xyCi6Xi5NOOomPPvoIgM2bN5OTk0NsbCznnXcexcXFfm3v3//+N2PHjsXlcjF16lS2bdvmm/ef//yH7Oxs4uPjueKKK6iqqupVrcHA71uQrVYrK1eu5LrrruM3v/kNZ5xxBn/6058AGD9+PI888ggXXXTRIetkZWXx3HPPcccdd1BWVsa5557LI488AsCJJ57IE088wcKFCzlw4ABXXnllj/2umKXj7cftMhIj2bijQrchi4iYraUBHkgx57Pv3APhnXfD0ZkFCxaQk5PDq6++ylNPPcXChQtZt24dM2fO5Cc/+QkrVqzggQce4MILL2TTpk2cddZZ5OXlHbGd3//+91x99dVcffXV3HbbbVx22WX83//9H4sWLeLvf/87xcXFnHfeefz2t7/lrLPO4uabb+a6665j1apVgfz2/c7vkAIwceJEPv/88yOmFxQUdLnOD3/4Q374wx92Ou/aa68N2rMn7QzD6DKkgM6kiIiI/xwOB83NzbhcLhYvXszixYt58cUXiYmJYfHixQD87ne/IykpiU2bNvHXv/6VpqamI7bTfretw+GgqamJ+Ph4nnzySd/8559/nqlTpzJv3jwAnnjiCUaMGEFpaSnDhg0bgG8aGL0KKcei8gPl1DbXYrPYSI9N90339ZVSrpAiImKqMJf3jIZZn90Ljz32GLfffjujRo1i9OjRPPDAAxQWFpKRkeFbxuFwkJKSwq5duzj11FO73d6KFSv4xS9+wdKlS8nOzua3v/0tkyZNorCwkMzMTN9yKSkpOBwOdu3aFVIhRQ8Y7EF7/yhp0Wk4bA7f9IzEtuf3VNTj8fTumqSIiASQxeK95GLG4Ec3HO08Hg9lZWW8/PLL7N+/n3nz5nH55ZczYsQIdu7c6VuusbGRPXv2MHLkyG6319DQQGtrK2vXrqWiooJp06Zx3XXXAZCWlsaOHTt8y+7evZumpqYetxlsFFJ60NmlHoDUeCd2q4XGFg97axrNKE1EREKI1Wpl7ty5LFu2jJKSEsAbXGbPnk1NTQ333nsvu3bt4qabbmL06NFMnjy52+21trYyc+ZMVqxYwb59+3zbA7jyyivZsGEDTz31FPn5+cyfP58LLrig225DgpFCSg98ISX+0JASZrMyMsF7mk/tUkRExB8vvPACL774ImPGjOHXv/41Tz/9NLGxsbz55pu89dZbjBs3jl27dvHqq69itXb/Ex0TE8OKFSu47777GD16NP/85z95/PHHARg5ciSrV69m+fLlTJgwAZfLxdNPPz0QXzGg1CalB+23H7f3kdJRRmIkO/fVk7+vnqlZiUfMFxER6Wj69OmdPvPupJNOYuPGjb3e3iWXXMIll1zS6bzvfe97nd7sEkp0JqUbHe/sGRM35oj5usNHRESk/yikdKOkvoSG1gbsVjsjY45sbJSRpJAiIiLSXxRSutF+FiU9Jp0wa9gR83UmRUREpP8opHSjqzt72mW23YZcuL+BFrdnwOoSERE5FiikdKO7RrMAw2IcOMNsuD0GRfsbBrI0ERGRQU8hpRvdNZoF71OddclHRESkfyikdMFjeNhZ5e0BsKszKaDGsyIiIv1FIaULu2t30+huJNwaTlp0WpfLZbadSdmpkCIiIhJQCild2Fa1DYDMuExsVluXy/ku9+hBgyIi0s8KCgqw9OJ5QaFOIaULPTWabac2KSIiIv1DIaULPd1+3K49pOytaaS+qbXf6xIRETlWKKR0wd+QEucKJyEyHICCCp1NEREZaIZh0NDSYMpgGEavas3Ly+OUU07B5XJx0kkn8dFHHwGwYcMG34MAJ0+ezJYtW3zrrF69mqysLIYMGcIzzzwTyF0X9PSAwU60elrJr84Her7cA96zKfvrm8nfV8+JKbH9XZ6IiHRwoPUAp/71VFM++6MrPsIV5vJ7+QULFpCTk8Orr77KU089xcKFC/n444+ZM2cOCxcuZN68eSxdupTbbruN119/ndLSUi677DKWLVvG9OnTufLKK/vx2wQfnUnpRGFtIS2eFpx2JyOiRvS4vBrPioiIPxwOB83NzbhcLhYvXswnn3wCwGeffcYdd9xBUVERtbW1fPPNNwC8+eabZGZmcsMNN5CVlcU999xjYvUDT2dSOtHeaDYzNhOrpeccp8azIiLmcdqdfHTFR6Z9dm889thj3H777YwaNYrRo0fzwAMPcO6557Js2TL++Mc/kpGRQVpaGm63G4CSkhLS0g52g5GZmRnQ+oOdQkon/G2P0k59pYiImMdisfTqkotZPB4PZWVlvPzyy9jtdv7whz9w+eWX8+qrr7J8+XK2b9/OsGHDeP311/n0008BGDp0KHv27PFto7Cw0KzyTaHLPZ3YXtm7kNLe6+zO8rpeN6ISEZFjg9VqZe7cuSxbtoySkhLAG1yqqqoAqK6uZsOGDdx6662+35KZM2fy9ddf8+yzz7Jjx45j7nKPQkon/O0jpV36EG9IqWlspbKhpd/qEhGR0PbCCy/w4osvMmbMGH7961/z9NNPc9555zFr1iwmTpzIjTfeyLx589izZw+lpaWkpqayYsUK7r33XnJycpg6darZX2FAWYwQ/qd/TU0NsbGxVFdXExMTE5Btuj1uznzpTKqbqnn7krdJjkz2a73Tl77DnupG/jF/CiePSghILSIicqTGxkby8/PJyMggIiLC7HKkC939nfz9/VablMPYrDbev+x99tbvZZhrmN/rZSRFsqe6kZ3l9QopIiIiAaDLPZ2wWCwMjxreq+cj6A4fERGRwFJICZCMxChAIUVERCRQFFICJFNnUkREBlQIN6k8JgTi76OQEiAdL/d4PPoPR0Skv9hsNgCam5tNrkS609DQAEBYWFift6GGswGSGu/EbrXQ1OqhpKaREXG964VQRET8Y7fbcblclJeXExYWhtWqf28HE8MwaGhooKysjLi4OF+o7AuFlACx26yMHOJiZ3k9+eX1CikiIv3EYrEwfPhw8vPz2bVrl9nlSBfi4uJITvavG4+uKKQEUGZipDek7KsjZ0yi2eWIiAxa4eHhjBkzRpd8glRYWNhRnUFpp5ASQBl6ho+IyICxWq3qzG2Q04W8ANJtyCIiIoGjkBJA6tBNREQkcBRSAiiz7WnIRfsbaG71mFyNiIhIaFNICaCh0Q5c4TY8BhTubzC7HBERkZCmkBJAFotFl3xEREQCRCElwA6GlDqTKxEREQltCikBpmf4iIiIBIZCSoBlJnlvQ95aUmtyJSIiIqFNISXATslIAODL3dXUNLaYXI2IiEjoUkgJsBFxTtKHuHB7DD7aud/sckREREKWQko/mJrlfW7Phu37TK5EREQkdCmk9AOFFBERkaOnkNIPpmQOwWKBbWV1lNU0ml2OiIhISFJI6QfxkeGcmBIDwIYdOpsiIiLSFwop/eTgJZ8KkysREREJTQop/SSnQ7sUwzBMrkZERCT0KKT0k0mjEgi3WSmpblTvsyIiIn2gkNJPnOE2Th4VD+guHxERkb5QSOlHU7OGAGqXIiIi0hcKKf2ovfHsxh37cHvULkVERKQ3FFL60bgRsUQ77NQ0tvLVnmqzyxEREQkpdrMLGMzsNiunjR7C21tKWb99H+NT48wu6djjcYO7Gdwt4Glte23pMN4+z+0dNzxguL3jvveew6a3zfN4Dr73DYb3FaPDuHHY+OHz28ahw7LtX6Dj+saR03zrQIeVepjWYfoRd54dNt7TnWlHdefaUZ5d1F1zIgPje3dCRIwpH62Q0s+mtoWUjdsrWDA9y+xygktLIxyobBv2w4EqaDkArQe8ry0Nba8dhwZobewwr8G7nZYD3umeFnC3HgwiR/tDKCJyrMu5RSFlsMoZ422XsqlgP40tbiLCbCZX1E9aGqGqEOrLDgaPhv0dQshhQ8N+bxgZcBawhYE1DGx2sIV731vtYLWBxdr2amt7tXR4bz3sfYdlLZaD07C0jVs6jFt7GG+rDbzT27fRXnN303zrdNhGj9PofhqdTevlst2u4+e6/anbukTEJ9xl2kcrpPSz0UlRDI12UFbbRO6uSk5va0wbcgwDGipgfz5UFkBl+2uBd1rtnr5t12IDZ3zbEAdhrrbBedjgAntE9/PsjrbQYe8QRNoCiC287f0gDYkiIoOQQko/s1gs5GQl8nLebtZv3xfcIcXjaQsf+R3CSIehua779cOjIHr4wdDhSugQQOI7n+6I0b9oRUSkUwopA+D0tpCyYUeQ9Zfi8UDZV1Cw/uDQWNXNChaISYH4dIjPaHtNh4S2964hChwiIhIwCikDoL1Tty+Lq6g+0EKsM8ycQvwJJfYIbwBpDx4dw0jcSAiLGPi6RUTkmKSQMgCGxzrJTIpkZ3k9H+6s4PsnJg/MBx8eSnZt8DZa7SgsEkZNgfQcSJ8Gw0/ytt0QERExWa86c8vNzSU7Oxun08nMmTMpKyvrcZ0NGzYwbtw4oqOjmTVr1iHrvPXWWxx33HFERkZy+eWX09DQ0PtvECI6PhW5X5V/Ax8+AS9eCQ9lwhM58Ob/B1+v9gaUsEjIOhvOvgdueAf+v11w1T+8t5ilTlJAERGRoOF3SPF4PMyZM4fZs2ezbds2nE4nt9xyS7frNDY28sMf/pDbb7+drVu3ArBo0SIAqqqquPTSS7n11lv56quvKCgo4P777z+KrxLcTh/djyHFMGDb2/Ds+fDoZHjz5wolIiIS8vy+3PPee++xf/9+7rnnHux2O4sXLyYnJ4f6+noiIyM7XSc/P59JkyZxzTXXADBr1iz+8pe/ALBq1SpSU1O58cYbAbjjjju45ZZbBm1QmZI5BKsFdpTXs7e6keTYALTtaG2GL1fCB8uhbIt3msUGmWd6L92kT4OUbAUREREJSX6HlA0bNjB58mTsdu8q2dnZuN1ucnNzmTZtWqfrnHDCCbz22msA1NbW8sorrzBhwgTf9k4//XTfsqeeeiqFhYUUFRWRlpbW5y8UrGJdYYwbEcvnxdVs2L6POSen9n1jB6rg02fgoyegtsQ7LTwKJv4ITpsPcYNv/4mIyLHH75Cyd+9eEhMP9vFhtVqJj4+ntLS0x3ULCwvJyMhg9OjR/O1vf/Ntb9y4cb5lhgzx3gFTWlraZUhpamqiqanJN15TU+Nv+UFhalbi0YWUqiL48HHIffZgnyVRyXDajXDydd7O0ERERAaJXjWcNQ57oJdhGFj86BcjJSWF9evXExsby1133dXp9trfd7e9pUuXEhsb6xtC7YzL1PbGszv2HbEvu1XyBfxjHiw7CT581BtQkk6ACx+Dm7/0ti9RQBERkUHG75AyfPhwysvLfeNut5uqqiqSk7u+nbaoqIitW7dit9uZMmUKd955JytWrOh0exUV3o7OutveokWLqK6u9g1FRUX+lh8UTh4Vj8NupbSmiR3lPfTeahiw/d/w3IXw5DT48iXvE3czzoAr/w4LPoAJV4I9fGCKFxERGWB+h5Rp06axadMmWltbAcjLy8Nut/vamHTmjTfe4Prrr/eNWywWbDabb3sbN270zfvwww9JT09nxIgRXW7P4XAQExNzyBBKIsJsTEqPB2DD9i56n3W3wGcvwONT4fk5sHOdtzHsdy+BH6+DH/0Txpyjnl1FRGTQ8zuk5OTkkJSUxOLFiykuLmbJkiVcfPHFuFwuampqaGlpOWKdGTNmkJeXx6pVq9i1axePPvoo5513HgAXXnghJSUlPP744xQUFPDQQw9x1VVXBe6bBan2Sz7rO7sVed92+MN0eOVGbydsYZFw2gL4nzy45E+Q0nUgFBERGWz8DilWq5WVK1eyevVqsrKyaGxs5OGHHwZg/PjxrFmz5oh1srKyeO6557jjjjsYP348CQkJPPLIIwDExsby0ksv8fDDDzN27FgyMjJ8fagMZlPb+kv5cGcFrW7PwRmbX/YGlNLN3mfgnHU33PoVnLsU4keZU6yIiIiJLEavWnAGl5qaGmJjY6murg6ZSz9uj8GEJf+iprGVVQtOZ0KKC/71C/j4D94FRk2FOX+CmOHmFioiItJP/P397tXdPXL0bFYLU0a3PXBw8xfw5+8fDCg5t8I1rymgiIiIoAcMmiInKxH31te5eNOTYNRBRBxc/Ac47vtmlyYiIhI0FFIGmruF80sf5+rwJ8AAT8rJWC99BuJGml2ZiIhIUNHlnoFUvRuemUXcZ08A8KfW/8fGM/6igCIiItIJhZSBsv0db6dsRR+BI4ZnUpfwf61Xsz6/1uzKREREgpIu9/Q3jxvW/RL+8xBgQPJ4uPRZYneFw/bP2dBZfykiIiKikNKvakvhH9dDwfve8ZOvg3N/CWERnG5vBGDznmqqGpqJc6l7exERkY50uae/5L/vvbxT8L6359iLn4LzH4GwCACGxUQwZmgUhgEf7Oiii3wREZFjmEJKoHk88P5v4LkLoK7U+7TiH6+F8ZcesWjHpyKLiIjIoRRSAqm5Af5+LbyzBAwPnDQX5r0DSd/pdHFfSOnqYYMiIiLHMLVJCZTavfDC5bAnD6xhMOvXMPFH3T6t+NTMBKwWyN9Xz+6qA4yIcw5gwSIiIsFNZ1ICoeQLeGqGN6A44+GaV+Hka7sNKAAxEWGclBYHoLt8REREDqOQcrS+XgN/PhdqdkPicTDvXUif6vfq7U9F3qiQIiIicgiFlL4yDNiwDF68ElrqIfN7cP3bkJDZq80cbDxbQQg/kFpERCTgFFL6orUZXvtvePtuwIBJ/wVXrgRnXK83NXFUHBFhVsprm9hWVhfwUkVEREKVQkpvNeyHv/wA8p4HixXOfRBmPQy2sD5tzmG3cUp6AgDrt+mSj4iISDuFlN7Ytw3+eBbsWg/h0XDFS3DajT02kO1J+yWfjeovRURExEe3IPtr5zp46RporIbYkXDF32DY2IBsOqctpHy4cz+tbg92m7KjiIiIfg398cnT8JeLvQEl7VTvHTwBCigAY4fHEOcKo66plc+LqwO2XRERkVCmkNIdjxvevBNW3wyGG8ZdCte8BlFJAf0Yq9XC6aOHAOovRUREpJ1CSleaauGFufDho97x7/0CLv6D7wGBgXb66PYu8hVSREREQG1SOldVCH+9HMq+AnsEXPQ4fPfifv3I9nYpuYWVNDS34grXn0ZERI5tOpNyuOrd8NRZ3oASNQyufb3fAwrAqCEuRsQ5aXEbbCqo7PfPExERCXYKKYeLSYGss2DYOG8D2dSTB+RjLRYLU7PULkVERKSdrikczmKB85eBuwUcUQP60VOzEnnpk2KFFBEREXQmpXN2x4AHFDjYeParPTXsr28e8M8XEREJJgopQSQp2sHxydEAvJK32+RqREREzKWQEmSunjIKgD/8ZydNrW6TqxERETGPQkqQueTkVIbHRrC3ppGVnxSbXY6IiIhpFFKCjMNu4ydnZALw+LodtLg9JlckIiJiDoWUIHT55JEkRjnYXXWAVblqmyIiIscmhZQgFBF28GzKo+u206qzKSIicgxSSAlSV542koTIcHZVNPDPL/aYXY6IiMiAU0gJUq5wO9fnZACw/N3tuD2GyRWJiIgMLIWUIHbNlFHEOsPYUV7PG5tLzC5HRERkQCmkBLHoiDCum5oOeM+meHQ2RUREjiEKKUHuutMziHLY+XpvLW9vLTW7HBERkQGjkBLkYl1h/Oh0by+0v393G4ahsykiInJsUEgJAdfnZOIKt7F5dw3rvik3uxwREZEBoZASAhIiw7nqNO/ZlN/pbIqIiBwjFFJCxA3TMnDYreQVVrFhe4XZ5YiIiPQ7hZQQMTQ6grmTRwLesykiIiKDnUJKCPnJmZmE26x8nL+fj3bqbIqIiAxuCikhZHisk0smpQLw+3e3m1yNiIhI/1JICTHzzxyN3Wph/fZ95BZWml2OiIhIv1FICTFpCS5+MGEEAL9/R21TRERk8FJICUELv5eF1QJrvynny+Jqs8sRERHpFwopISg9MZILTkoBvL3QioiIDEYKKSHqv2dkYbHAv7aUsrWkxuxyREREAk4hJURlDY3mvO8OB2D5Wt3pIyIig49CSgj77xlZALz+ZQnby2pNrkZERCSwFFJC2AnDYzhn7DAMAx5du8PsckRERAJKISXE/c+MMQC8+tlu8vfVm1yNiIhI4CikhLhxqbFM/04SHgMeU9sUEREZRBRSBoGftp1NWZW3m6L9DSZXIyIiEhgKKYPAyaPimZo1hFaPwePvqW2KiIgMDgopg0T72ZS/f1JMSfUBk6sRERE5egopg8RpmUOYnJ5As9vDcj0hWUREBgGFlEHkprO9Z1NWfFTI8x/uMrkaERGRo6OQMohMzUrkp20dvN316mZe/7LE5IpERET6TiFlkLn1nOO44tSRGAbc/OJnbNi+z+ySRERE+kQhZZCxWCz834Xf5f99N5lmt4cfP/cJXxRXmV2WiIhIr/UqpOTm5pKdnY3T6WTmzJmUlZX1uM6nn37KxIkTiYmJ4cILL2T//v2+ecnJyVgsFt8wadKk3n8DOYLNauGRy7M5ffQQ6pvdXPv0JnaW15ldloiISK/4HVI8Hg9z5sxh9uzZbNu2DafTyS233NLjOldccQUzZ87kiy++oKysjLvvvts3v6qqiq1bt1JZWUllZSXr1q3r8xeRQznsNv5wzSTGjYhlf30zV//pY/ZWN5pdloiIiN8shmEY/iy4du1aLrroIioqKrDb7eTm5pKTk0N5eTmRkZGdrrN9+3bGjBlDXV0dkZGRPPbYYzz55JN8/vnnNDU1ERERQVNTE+Hh4X0qvqamhtjYWKqrq4mJienTNga7fXVNXPrEB+zcV89xw6J46SdTiHP1bX+LiIgEgr+/336fSdmwYQOTJ0/GbrcDkJ2djdvtJjc3t8t1XC4XjzzyiC/EVFRU4HQ6Ae9ZlIiICObOnYvT6eSMM85g9+7d3dbQ1NRETU3NIYN0LzHKwbP/NZlhMQ6+La3jv57ZxIFmt9lliYiI9MjvkLJ3714SExMPrmi1Eh8fT2lpaZfrpKSkcNNNNwFQXV3NH//4R66++mrAG1IaGxs566yz2LJlC1arldtuu63bGpYuXUpsbKxvSEtL87f8Y1pagovn/utUYiLs5BZWMX/Fp7S4PWaXJSIi0q1eNZw9/MqQYRhYLJYe16uvr+f8889nwoQJzJ8/H4DMzEz27NnDggULyMjIYP78+axdu7bb7SxatIjq6mrfUFRU1Jvyj2nfSY7m6etOISLMyrpvyrnj71/g8fh1pU9ERMQUfoeU4cOHU15e7ht3u91UVVWRnJzc7XoNDQ18//vfx+Vy8eKLL2K1ej8yLCyM4cOH+5aLj4/v8fKNw+EgJibmkEH8d/KoBB6/8mRsVgur8nZz35qtRwRPERGRYOF3SJk2bRqbNm2itbUVgLy8POx2OxMmTOh2vQULFhAZGcmrr75KRESEb/qjjz7KjBkzfOOFhYWkp6f3snzpre8dP5SHLhkPwJ835PPYOj01WUREgpPfISUnJ4ekpCQWL15McXExS5Ys4eKLL8blclFTU0NLS8sR63z++ee88sorLF++nAMHDlBVVUVVVRUA06dPZ+PGjaxevZpvvvmGZcuWce211wbqe0k3Lp6Yyi9mnQDAQ299w4sfF5pckYiIyJH8DilWq5WVK1eyevVqsrKyaGxs5OGHHwZg/PjxrFmz5oh1XnnlFaqrqznuuOOIj4/3DQAnnngiTzzxBAsXLmTatGmcffbZPfa7IoFzw7RMFkwfDcCdq77kzc17Ta5IRETkUH73kxKM1E/K0TEMg0Uvf8mLm4oIt1t59rrJTBk9xOyyRERkkAt4Pyky+FgsFu676Lt8/8RhNLd6mPfcJ2zeXW12WSIiIoBCyjHPbrOy7PIJnJqRQF1TK9c+/THvbyvveUUREZF+ppAiRITZeOpHkzgxJYZ9dd7n/Mx77hN2VdSbXZqIiBzDFFIEgJiIMF748WlcNzUdm9XC21tKOefh//Dgm19T19RqdnkiInIMUsNZOcK20lqWrN7C+9v2ATA02sHPzz2eH0wYgdXacw/DIiIi3fH391shRTplGAb/3lrGfWu2sKuiAYDstDgWnz+WCSPjTa5ORERCmUKKBERTq5s/ry9g+bvbqG97evKcian8/NzvMDQmooe1RUREjqSQIgFVVtPIg29+wz9yiwGIDLfx3zPG8F856TjsNpOrExGRUKKQIv0ir7CSe/65hc+LqgAYNcTFL2aN5ewThvr1RGwRERGFFOk3Ho/Bqrzd/PLNrymvbQJg2phE7p49ljHDok2uTkREgp1CivS7uqZWHl27nT+9n0+z24PNauHCk1L4/neTOWNMEs5wXQYSEZEjKaTIgNlVUc99a7by9pZS3zSH3cq0MYmcM3YYM44fRlK0w8QKRUQkmCikyID7pGA/q78o4e0tpeyuOuCbbrHAxJHxnH3CMM4ZO4ysoVEmVikiImZTSBHTGIbB13treXtLKW9vKeXLwx5amJkYyTljvYFlwsh4bOogTkTkmKKQIkGjpPoA/95axttbSvlgxz5a3AcPuSGR4cw4fijnjB3GNLVjERE5JiikSFCqbWzhvW/LeXtLKWu/LqOm8eBzgcJtVkYPjWJM+zAsiqyh0Ywa4iLMpsdMiYgMFgopEvRa3B425e/nX22XhTq2Y+kozGYhfUikL7SMGRpF1tAoMhIjiQjTmRcRkVCjkCIhxTAMCvc3sK20jm1ldWwrq2VHmfd9Q1t3/IezWmDUkEiy2s68jExwMSTKwZCocBIjHSRGh+MKtw/wNxERkZ4opMig4PEYlNQ0sq20lu1ldWxvCy7bSmsPuVTUFWeYjSFR4QyJcpAYGX7wfZSDxKhwhkR6Q82QyHCiIuw4w2zqOVdEpJ8ppMigZhgG5bVNB0NLWS0lVY3sq2+moq6JfXVNNLZ4er1di8UbbFzhdiId3ldXuA1XuI3I9veO9vcHx51hNsJsVsJsVhx2a9t7C2F2K+E2K+Edph067p2mYCQixxJ/f791LlxCksViYWhMBENjIjg9K7HTZeqbWqmoa2ZffRMVdd7wUlHfTHmt97Wirm16fRP765vxGGAY0NDspqHZzb66gfs+NqsFm8XifbVasFrAbrNitViwWcFutWK1csgyNqsVW9s0i8WCxQJWi3ddi8WChbZxq/cVDp3fcTmLBbzv2t63j3v/d9hytC3XtkbbhI7rd5h8xLY7vh5c6vDpnc3tbH734a677He0sfBYDZbH6Nc+pt16znFER4SZ8tkKKTJoRTrsRDrsjBzi6nFZj8egsdVNfZObhubWtqDS6huvb3LT0OKmoamV+uaDrweava+NLW6aWz20uD20uA1a3B6aWz00uw9Oax9vbj3yDI/bY+DGgM6b34iImGb+9NEKKSJmslotbZdv7ED/duFvGAZuj3FIcHF7DNyGgcdj0Orxzve0LecbOox72sZb294bBngMA6Nt+5728R5ePQYYeMcNb3Ft2/BuBzg4TodpHdZrX6Z9unfc6PB9j/z+nU4/Yj8dPt/odn5n2+hqwaO9xm3mRfLD94NIfzPzBgSFFJEBZrFYsNss2G2o8zoRkW6ohywREREJSgopIiIiEpQUUkRERCQoKaSIiIhIUFJIERERkaCkkCIiIiJBSSFFREREgpJCioiIiAQlhRQREREJSgopIiIiEpQUUkRERCQoKaSIiIhIUFJIERERkaAU0k9Bbn/ke01NjcmViIiIiL/af7fbf8e7EtIhpba2FoC0tDSTKxEREZHeqq2tJTY2tsv5FqOnGBPEPB4Pe/bsITo6GovFErDt1tTUkJaWRlFRETExMQHb7mCn/dY32m+9p33WN9pvfaP91jfd7TfDMKitrSUlJQWrteuWJyF9JsVqtZKamtpv24+JidEB2Qfab32j/dZ72md9o/3WN9pvfdPVfuvuDEo7NZwVERGRoKSQIiIiIkFJIaUTDoeDxYsX43A4zC4lpGi/9Y32W+9pn/WN9lvfaL/1TSD2W0g3nBUREZHBS2dSREREJCgppIiIiEhQUkgRERGRoKSQcpjc3Fyys7NxOp3MnDmTsrIys0sKCaeddhoWi8U3JCYmml1S0CopKeHMM8/ks88+803Tcdezzvabjrvu7dy5kzPPPJPo6GimT5/Orl27AB1v3elqn+lY697WrVuZMmUKUVFRTJs2jW3btgFHf6wppHTg8XiYM2cOs2fPZtu2bTidTm655RazywoJVVVVvPXWW1RWVlJZWcnOnTvNLiko/eQnPyElJYX//Oc/vmk67nrW2X4DHXc9+fGPf8zIkSPZvHkzQ4YMYeHChTreetDZPgMdaz254ooruOCCC/j22285/vjjufHGGwNzrBni8+677xoxMTFGS0uLYRiG8emnnxpOp9Ooq6szubLgl5ycbHz77bdmlxH0ysvLjfz8fAMw8vLyDMPQceePzvabYei4605TU5NhsViMr776yjAMw1izZo0RExOj460bXe0zw9Cx1p39+/cbOTk5RlNTk2EY3v2WnJwckGNNZ1I62LBhA5MnT8Zu9z4tIDs7G7fbTW5ursmVBb+qqiruvvtunE4n2dnZfPXVV2aXFJQSExNJT08/ZJqOu551tt9Ax113Wlpa+NWvfkVGRgYAFRUVOJ1OHW/d6GqfgY617sTHx/P+++8THh5Oc3MzL730EhMmTAjIsaaQ0sHevXsPuc5otVqJj4+ntLTUxKqCX3NzM42NjWRlZbFlyxbGjh3LDTfcYHZZIUPHXd/ouOteZGQkt912G06nk5aWFn73u99x9dVX63jrRlf7TMea/1wuF2+88QbLly8PyLEW0g8Y7A/GYX3bGYYR0CcsD0ZhYWEUFxczYsQIAG666SZOO+00Dhw44PtXiHRPx13v6bjzT2trK1deeSVWq5UlS5Zw++2363jrweH7TMea/z7++GP+93//l//5n/8hPT39qI81nUnpYPjw4ZSXl/vG3W43VVVVJCcnm1hV8LNYLL7/eMF76g+gtrbWrJJCio67vtFx1zOPx8Pll1/O9u3beeONN3A6nTreetDZPtOx1r3y8nLy8vIAmDhxIg888ABr1qwhKSnpqI81hZQOpk2bxqZNm2htbQUgLy8Pu93OhAkTTK4suK1evZrRo0f7xgsLC3G5XCQlJZlYVejQcdc3Ou56tmTJErZv3867775LQkICoOOtJ53tMx1r3cvLy2PWrFm+8fYzJdOnTz/6Yy2gTXxDnNvtNjIzM40777zTKCoqMs4//3zjqquuMrusoFdWVmZERUUZTz31lLFz505jxowZxvz5880uK6jR4S4VHXf+67jfdNx1r6SkxIiJiTE2btxoVFZW+gYdb13rap/t3btXx1o3Kisrjfj4eGP58uVGUVGRcfXVVxvTpk0LyLGmkHKYTz/91Bg/frzhcDiMc845xygrKzO7pJCwZs0a47jjjjPi4uKMq6++2qitrTW7pKDGYbfS6rjzz+H7Tcdd15555hkDOGLIz8/X8daF7vaZjrXurV271hg/frwRFRVlzJw50ygoKDAM4+j/v01PQRYREZGgpDYpIiIiEpQUUkRERCQoKaSIiIhIUFJIERERkaCkkCIiIiJBSSFFREREgpJCioiIiAQlhRQREREJSgopIiIiEpQUUkRERCQoKaSIiIhIUPr/AaS2odjVWaNCAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import random\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "\n",
    "transfer_matrix = np.array([[0.6, 0.2, 0.2], [0.3, 0.4, 0.3], [0, 0.3, 0.7]],\n",
    "                           dtype='float32')\n",
    "start_matrix = np.array([[0.5, 0.3, 0.2]], dtype='float32')\n",
    "\n",
    "value1 = []\n",
    "value2 = []\n",
    "value3 = []\n",
    "for i in range(30):\n",
    "    start_matrix = np.dot(start_matrix, transfer_matrix)\n",
    "    value1.append(start_matrix[0][0])\n",
    "    value2.append(start_matrix[0][1])\n",
    "    value3.append(start_matrix[0][2])\n",
    "print(start_matrix)\n",
    "\n",
    "#进行可视化\n",
    "x = np.arange(30)\n",
    "plt.plot(x,value1,label='cheerful')\n",
    "plt.plot(x,value2,label='so-so')\n",
    "plt.plot(x,value3,label='sad')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffef4a01",
   "metadata": {},
   "source": [
    "可以发现，从10轮左右开始，我们的状态概率分布就不变了，一直保持在[0.23076934,0.30769244,0.4615386]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9678a82",
   "metadata": {},
   "source": [
    "# Metropolis-Hastings采样python实现  \n",
    "https://zhuanlan.zhihu.com/p/37121528\n",
    "\n",
    "假设目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵 $Q(i,j)$ 的条件转移概率是以 $i$ 为均值,方差1的正态分布在位置 $j$ 的值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "92755ed2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGbCAYAAADawqrfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABjAUlEQVR4nO3de3jT5fk/8PcnTQ9JD+kR2tJCCqUoICsHYWKLCApOYWzgHDrYjzk3Obg5+DImm9qCOpxsDDZRt/Hd0KFfJ27AKAgeoAVaUaTljFAOpUCP0Calh6Rtkt8ftZHS0+dJ88mp79d15bqa5HmSO6Ekd5/D/Ug2m80GIiIiIi+icncARERERKKYwBAREZHXYQJDREREXocJDBEREXkdJjBERETkdZjAEBERkddhAkNERERehwkMEREReR21uwNQitVqRUlJCUJDQyFJkrvDISIiIhlsNhtu3LiB+Ph4qFSdj7P4bAJTUlKCxMREd4dBREREDrh8+TISEhI6vd9nE5jQ0FAALW9AWFiYm6MhIiIiOWpqapCYmGj/Hu+MzyYwrdNGYWFhTGCIiIi8THfLP7iIl4iIiLwOExgiIiLyOkxgiIiIyOv47BoYIiJvYrPZ0NzcDIvF4u5QiBTl5+cHtVrd4xInTGCIiNyssbERpaWlqK+vd3coRC6h1WoRFxeHgIAAhx+DCQwRkRtZrVZcvHgRfn5+iI+PR0BAAItvks+y2WxobGxEZWUlLl68iMGDB3dZrK4rTGCIiNyosbERVqsViYmJ0Gq17g6HSHEajQb+/v64dOkSGhsbERQU5NDjcBEvEZEHcPSvUCJv5Izfd/6PISIiIq8jlMDk5+cjNTUVGo0GU6ZMQUVFRbd9Dh8+jFGjRiEsLAwzZsxAVVWV/b7du3cjJSUFwcHBmD17dpsFbF3dR0S906krNdA/s6PTS96X19wdIhG5iOwExmq1YtasWZg2bRoKCwuh0WiwePHibvs89thjmDJlCo4dO4aKigo8//zzAACDwYBHHnkES5YswcmTJ1FUVISXXnqp2/uIqHc5cKrSnqA8+Or+Lts+tvEze9tNOYUuirD32bhxIyRJanfJzs52d2iytb4GPz8/9O/fH8uWLUNjY2ObNtnZ2dDr9S6LqaPn27hxIyZOnKjI8xUVFXn1gnHZCUxOTg6qqqqQmZmJhIQEZGRkYMuWLairq+u0z4ULF3D27Fk899xz0Ov1mDt3Lvbvb/kA2rJlCxISEjB//nzo9XosW7YMmzZt6vY+IuodjhQZoH9mB+a89blD/Z/94Cz0z+zAfw8WOzkyz2Sx2vDp+evYduQqPj1/HRarTbHneuyxx1BdXY19+/YBAKqrq1FdXY20tDTFnrM7EydOxMaNG4X6DB8+HCUlJfjzn/+Mt99+GwsXLmxzf1paGo4dOyYci16vdyiZc/T5utLV+9K/f39UV1c79flcSfYupNzcXIwdOxZqdUuX1NRUWCwW5OfnIz09vcM+Wq0Wa9euRXBwMADg+vXr0Gg09scbP368ve24ceNQXFyMy5cvd3lfYmJih89lNpthNpvt12tqauS+NCLyIFerGnD3K3uc9ng/33ocP996HPuW3ov+0b65y2fXiVKs2H4KpUaT/bY4XRAypg/FA8PjnP58AQEBCAgIsJ8WHB4e7vTncAU/Pz/07dsXM2bMgEajwUMPPYRXXnkFkZGRAAC1Wu3Sw4Bd/Xwqlcpr/+0AgRGYsrIyREdHf91RpUJERATKy8s77RMfH4+nn34aAGA0GrFhwwbMnTu3w8eLiooCAJSXl3d5X2dWrVoFnU5nv3SW6BCR5xr0zA6nJi83m/D7vdA/s0ORx3anXSdKsWBTfpvkBQDKjCYs2JSPXSdKXR5Tbm4uRo4cCa1Wi7Fjx+LUqVP2+zIzMzFv3jzs2rULqampWLZsmf2+4uJijB8/HjqdDvPnz8eYMWPwy1/+EgBw7tw53HfffQgLC8Po0aNx6NAhAMD8+fMhSRJycnLwox/9CJIkYf78+cIxT5o0CZIkoaCgwH5bZ1NIH3/8MYYOHQqtVou7774bhYUt05UPPPAAJEnCpUuXcO+990KSJLz88sv2fvPmzUNmZiY2bdqEIUOG4NVXX23zuJ09X2NjI2bOnImQkBA8+OCD9vWnHbWXJAlFRUWy3pfOppBOnDiBtLQ06HQ6PPjgg7hy5Uqb5/vvf/+LAQMGIDIyst1rcCWhRbw2m63ddTnzZ3V1dZg+fTpGjhyJBQsWdPh4rT+3Pl5X93Vk+fLlMBqN9svly5dlvCIi8hT6Z3bAFUX0fSmJsVhtWLH9FDqaLGq9bcX2U4pOJ92qdb3kzJkzceHCBdx1111YunRpmzYnTpzAsmXL8Nxzz+HJJ5+03/7MM89g5MiROHz4MHbt2oVnnnkGv/rVr9Dc3Ixvf/vbmDt3Lk6cOIG5c+fi0Ucfhc1mwx//+EdUV1fj7rvvxvr161FdXY0//vGPwnGr1WpER0ejsrKy27Zz587Fj3/8Y5w9exbDhw/H8uXLAQD//ve/UV1djcTERGzfvh3V1dXt1oru3r0br732GtasWYPvfOc7smL79NNPMXr0aBw/fhxqtRpPPfVUt30cfV9qa2sxZcoU3H///Th27BgSExMxY8YMWK1WAC0zKb/73e+wc+dOrFy5Ev/zP/8Dk8nUzaMqQ/YUUlxcHE6fPm2/brFYYDAYEBsb22W/+vp6TJ06FSEhIXj33Xfte7/j4uLa/KJcv34dABAbG9vlfZ0JDAxEYGCg3JdDRB7E1UmF/pkdKHr5IZc+pxI+v1jVbuTlZjYApUYTPr9YhbsGRbksriNHjiAiIgLHjh3DjRs3cObMmTb3nzhxAl9++WW70YMjR45g7dq1SE5Oxl133YXCwkI8/PDDyM3NxZdffmkf0QdaRvVLS0sRHx8PjUYDtVoNrVbboykRSZLa/aHekcDAQJjNZkREROAvf/mL/fbW5RIqlQohISEdxtK6NlSn08mOKz4+Hr/+9a8hSRIyMjLwzW9+s9szszQajUPvS1ZWFsLCwpCRkQEA+NOf/oSYmBj7iFdtbS3eeOMNDBs2DIMHD8bPfvYzlJeXY8CAAbJfj7PIHoFJT0/HoUOH0NzcDAAoKCiAWq3GyJEju+y3cOFCBAcHY9u2bW2q7aWnpyMvL89+/eDBg9Dr9ejXr1+X9xGRb3HXiIgvjMRU3JD3l6/cds6gUqmwbt06JCQkYNGiRTAaje2+bKdPn97hVElycjI+/fRTmEwmHD16FEOHDgUAXLlyBQMGDMCRI0fsl4sXLyImJsZpcVssFly7dg19+/bttu3bb7+Njz76CLGxsUhPT8cXX3wh+3l++MMfCiUvADBgwAD7DERiYiKam5tx7Vr7kgFdbaqRq7i4GElJSfbrgYGBiI+Px6VLlwAAERERuOOOOwDAfo6RnKRPCbITmLS0NMTExCAjIwNXrlzBypUrMXPmTGi1WtTU1KCpqaldn6NHj2Lr1q149dVX0dDQAIPBAIPBAACYMWMGSktL8frrr6OoqAirV6/GnDlzur2PiHyHu5MIdz9/T/UJlVeCXW47Z8jJycGrr76KEydO4PPPP8ePf/zjdm1aRypuNWTIEPzxj39EaGgohg0bhunTpwMAEhISUFFRgT59+kCv1yMhIQGvvPKKfXQeaEmcevJFmpOTA0mSMHr06C7b1dfXo7m5GXv37sX169eRnp6OH/3oR23adBVLZ6+9K5cvX7Y/3tWrV+Hn54eoqChIktQmOewokRJ9XxITE3HhwgX7dZPJhJKSEvTv3x8AXLrIuDuyExiVSoXNmzcjKysLycnJMJlMWLNmDQBgxIgR2LGj/QfB1q1bYTQakZKSgoiICPsFAHQ6Hd577z2sWbMGQ4cORVJSkn0esav7iMg3eEry4ClxOGJsUiTidEHobHWghJbdSGOTIl0WU+sOUKPRiNzcXCxZskTWF6jZbMb69evxySef4OzZs3jvvffsSw7Gjh2Lfv36YenSpbh8+TJWrVqFHTt2tBmBGTRoEPbs2YPS0lJ8/PHH3U6xAC2jLuXl5cjKysK8efOwaNGibkdHmpubMWXKFLz99tv2UZDW9SE3x/Lhhx+itLQUn3zySbdxdOfKlStYvXo1ioqKsHLlSjz00ENQq9Xo168fSktLcf78edTV1dmnfW6NReR9mTZtGmpqarBixQpcunQJTz/9NAYNGoSxY8f2+HU4m9Ai3lGjRuHo0aMwmUz48MMP7b88RUVFHS5GysjIgM1ma3dpNWXKFBQWFqK+vh7vvvtum4PMurqPiLybpyUNnhaPXH4qCRnTW6ZZbk1iWq9nTB8KP5XripU98MADeOihhzBq1CjMnz8fP/nJT1BSUtLlLlKgZari/vvvx9SpU5GSkoKAgADMnDkTZrMZ/v7+2L59O7788ksMHToUWVlZ2LZtG/z8/Oz9n332WZw/fx5JSUn46U9/2i6p6MiJEycQHx+Pn//855g/fz5Wr17dbZ+wsDC8/fbbePHFFzFo0CBs374dr7/+eps2q1evxo4dO6DX6+3FW3vizjvvxIEDBzBixAjU1dXhtddeA9Ay5faLX/wCaWlpSEtLa7dYGhB/X0JDQ7Fr1y7s3r0bd9xxBy5duoRt27Z55Fldks1dk1cKq6mpgU6ng9Fo9KghL6LezpnJwvaFadh06Dz+dcg5W4XdsbDXZDLh4sWLSEpKcvhUXlfXgVHCnj178Itf/ALbtm2DTqfDhQsX8K1vfQt79uyxr7kg39HV773c72/Zu5CIiHrKGcnLicypCAn6+qPrd/1H4XezWn4+cKrS4cq9gPfuTnpgeBzuHxqLzy9WoeKGCX1CW6aNXDny0lPf+MY3MGjQIIwePRq1tbWIi4vDwoULMXz4cHeHRh6KIzBE5BLL3i/Ae1+UONz/41/cg+TYEFltRRKlDe+vaHNdAjD59lt2omzfLvvxRDljBIbI2zhjBMbzJrWIyOc0Nlt7lLwUvfyQ7OSltb3awU83G4Avy3gUCZGnYwJDRIpLefYDh/s6OqVz7rcPIXfZJIf6XqlugNU3B6eJfAYTGCJS1IJ/5nXfqBM9XY/SL1KD87990KG+e76s6NFzE5GymMAQkWIam6344GS1Q32dtZjWTyVh9cMjHOr7+YX21U6JyDMwgSEixTg6deTsnUDfG5OIYH/xHTk1ZotLD0IkIvm4jZqIFDHIwS3TSm1jPvnCgw5t4957pgL3KRCPLF+V0ncJwZ1Wn332GRYuXIizZ89i3LhxePPNNxU/r27jxo3YuHEjsrOzFX2eiRMnIicnB0FBQRg+fDief/55+7EGrebNmwe9Xo/MzExFY+nq+SZOnIh58+Zh3rx5Tn++zMxMFBUVYePGjU5/bGfhCAwROV2ZwYTuC7m3dyJzqtNjuZmjyVGZwXWHIXqDhoYGzJgxAwsWLMCpU6cQGhqKp556yt1hOdVvf/tbHDt2DA888AC++93vYt++fW3uf+211/DMM88IPWZ2dnaHh1jK4cjzdaWoqMh+QGRHnnnmGXvFX0/FBIaInO6bL4uf/5LSV9umQJ1Sdj6VLtzHkdfjy06dOoXq6mo88cQTSExMREZGRpdfht5Io9Fg8ODBeOGFF/DII49g3bp1be7XarUurdvj6ucLCgry+CN8mMAQkVNdrKhzqN+Hi+91ciQdG5rgWGHLe17+0MmReK/ExERIkoQXX3wRzc3NSE1NxX/+8x/7/du2bcOQIUMQHByMyZMno6SkpQaQXq/H/PnzERsbi1/96leYMWMGYmJicPjwYWRmZuLBBx/EpEmTEB4ejtmzZ9sPhuzOW2+9hcGDByM8PBw//OEPUV9fb7/vj3/8I+Lj4xEaGorZs2fDZBIfTZs6dSo+/7xthed58+a1mz6y2WxYtmwZYmJiEBERgaeeego2mw1lZWWQJAn33nsvLl26BEmSIEkSysrK7H0lScLJkyfx5JNPIjIyEkajsdvnA4DDhw9j0KBB6NOnD1566aVO22/cuBETJ04E0JKcJCUl2Z9XkiQcPHiwzeNmZmZ2ODW1fv166PV6xMfHIzMz03620rx58/Dcc89h0aJFCAkJwbBhw3DmzJkO309nYQJDRE5175ps4T6nVz7g/EC64MhU0iVDE2pNzQpE43369OmDTZs24fe//z2Sk5Pxz3/+035fVVUVvv/972P58uU4d+4cwsPD8eKLL9rvv3HjBlauXIlXXnkFjz/+OEaMGIFdu3YBAHbt2oXHH38cX3zxBS5duoTnnnuu21j279+P3/3ud3jnnXfw2WefoaysDL///e8BAGfOnMEvf/lLvPfeeygoKMD58+exYcMG4dcbFxeHysrKbtt9+OGH2LBhA/bs2YMDBw5g+/bt2LFjB/r27Yvq6mps374diYmJqK6uRnV1Nfr2bVvx+YknnkBYWBi2bNmC4OBgWbG9++672LRpE/79739jzZo12LlzZ7d9ysvLcfToUQCwx3LnnXd22+/f//43VqxYgY0bNyIrKwtvv/02/vznP9vv/8tf/oKwsDCcOHECffv2bZNQKYEJDBE5zf1r9gr3mTAoEpoAv+4bOtnR56cI9xmeuVuBSLzTww8/jEuXLmHevHn46U9/il/+8pcAWk4zvnTpEmbPno3z58+jsbGxzV/ic+fORUpKCvr27YsZM2YgMTERTU1NAIC7774bc+bMQXJyMn71q19h27Zt3cbx1ltv4fz587j//vsxbtw4HDhwAAcOHAAABAQEAADMZjP0ej0OHTrk0FodSZIg59SdwMBAWK1WNDY24vbbb8elS5cwbdo0SJKE8PBwhISEQKVSITw8HOHh4e2m3UaMGIHVq1fjnnvugVotbzr1iSeewF133YX09HQ89thjeP/997vto9Pp7CX6W2O5+WTvzvztb3/D4sWLMXHiRIwaNQorVqzAG2+8Yb8/MTERq1atgl6vx+zZs3H58mVZr8FRTGCIyClqTc0orKjvvuEt3vrJXQpE0z2d1h99QsTX3BRfE3+NvqakpATnz5+HTqdDZmYmPvjgA/zhD39AcXExAOA3v/kN+vXrh2eeeQbNzc2wWL5e0t26jqOj9RyJiYn2n+Pj41FeXt5tLFeuXMHPfvYzHDlyBEeOHMGpU6fw5ptvAgCSkpLwxhtvYOnSpYiKisIjjzyCigrxAoXl5eXtRks6MnHiRCxbtgxz5sxBTEwMfvKTn6CuTv6U6s9//nPh2FqngoCW96+0tOOT2UXi6ExxcTEGDhxovz5w4EBcunTJfv2ee+6x/xwQECAr6esJJjBE5BSjV4iPTuQ/e78Ckcj3+bPiu54m/F58lMnX/Otf/8ITTzxhvz5hwgSo1WpUV1fjnXfewZ49e3Dp0iUcOHCg3fbjrly8eNH+c3FxMeLi4rrtk5CQgGvXrkGv10Ov1+Pq1at46623ALQkWnfeeScKCgpQVFSE6upqvPDCCwKvtMWHH36Iu+++u9t2Fy5cwKxZs3D69GkcP34cn332WZsRCpVK1eWXutxpo5u1Jo0AcPXqVXuiJUlSm8Txiy++aNNPpWr5+hdJMhITE3H+/Hn79fPnz6N///72664+OJkJDBH1WK2pGWbBP7ZCA4DIkABlAhIwLilKuM/hC45VF/YVkydPRl5eHv7v//4PV69eRWZmJmJjY3HbbbfZF95WVVXhgw8+wAsvvCD7S/LgwYN48803UVhYiFdeeQUzZ87sts+cOXPwr3/9C//9739x/vx5LFmyxB7D8ePHMWXKFOTm5tpHIFoXnXanoaEB586dw8qVK7F582YsXbq02z4fffQRvvvd76KgoABms7nd8w0cOBAlJSUoKCjAuXPnUFBQICuWrvztb3/DZ599hry8PLzzzjuYMWMGAKBfv37Iy8uDzWbDp59+2maRNdCyrker1SIrKwuXLl1qt4i3I0888QTWrl2LnJwcFBQUIDMzE08++WSPX4OjmMAQUY85sjbk+EplCtaJCnVg6/asvzp+vpMvGDFiBP7xj38gIyMDKSkp2LNnD7Zu3YrAwEDMnTsXgwYNwu23344VK1bgySefxOnTp2Xt/pk+fTr+/ve/Y9SoURg0aBAyMjK67XPPPfdg/fr1WLJkCUaPHo077rjDvvh36tSpePLJJ/G9730PKSkpsNls+M1vfiPrNf7617/GHXfcgY8++ggffvghRo8e3W2fxx9/HBMmTMDUqVMxYsQIDB48GAsWLLDfn5CQgFWrVuG+++5DamoqcnNzZcXSlenTp+P//b//Z6/L893vfhcAsHDhQtTW1mLQoEF44YUX2r1uf39//PWvf8WTTz6J2267DVu2bOn2ub73ve/h2Wefxdy5c/Gtb30Ls2fPdmjay1kkm9KTVG5SU1MDnU4Ho9Ho8mEtot6kzGASrpOSu2wS+kVqFIpI0PTpaGy2Yl9h57tMnni4/Rfp7p9PwJD40B4/vclkwsWLF5GUlOTSOh+exhsqv5LzdPV7L/f7myMwRNQjjhR585jk5SsBahX8BAuxTf3Tvu4bEZFieBYSETmsssbc5vqG91d020d/YI9S4fTIPUNisOdLsR0q58pqkRwbolBEvYurzhQi38ERGCJy2J2//Vi4j6d+4askCWGC62HuW5ujUDRE1B0mMETkEGN9k3CfuwdFKxCJ84zRRwr3cfToBCLqGSYwROSQO18UPxvIHRV3RagkCTrBURhHjk7oiI/upyDqkDN+35nAEJGwWlMzGuWV07BLS45RJhgnG+3AKEyZQfyAwFb+/v4A0OYAQiJf1/r73vr77wgu4iUiYd9woO5LkL93/L2kkiREatWoqpd/cOM3X/7EoQMiAcDPzw/h4eH2EvdarbbdGTlEvsJms6G+vh4VFRWyz2DqDBMYIhJirG+Cpftmbdw7pI8isShlZP8ofPJl9+fw3Kz4Wj36R2sder7Y2FgAcOicHiJvFB4ebv+9dxQTGCIS8p114juP/FTeNaIgSUCf0EBU3DB33/grE36/1+FRGEmSEBcXhz59+thPZibyVf7+/j0aeWnFBIaIhFw0ii1+GevAmhJPMLyfTrguTEOjpUcLlf38/JzywU7UG3jHpDQReQRHtgyHaRxfpOdOKkmCv5/YyNHtz+9SKBoiuhUTGCKSTXTL8ITB3rHzqDPjHahbc7WqQYFIiOhWnEIiIlmqahuF+wSovftvJH8/FeSOwbQeo3Dm/RXod3vfjhtt3+6cwIiIIzBEJM+YFz8Sau8tdV+6MyFF7HUIlschIgcJJTD5+flITU2FRqPBlClTZG/5O3fuHIYOHQqDwWC/LTs7G5Iktbno9Xr7/bGxsW3uGzNmjEioROREVbWNwl/Mitd9mT69+4sT+PuJv45mK6vqEilN9v9Mq9WKWbNmYdq0aSgsLIRGo8HixYu77Td16lQMHjwYp0+fbnN7Wloaqqur7ZcXX3wRcXFx9vsNBgNOnz5tvz87O1v+qyIip5r+52yh9iEBvjW4u2/pvULtv7h4XaFIiKiV7DUwOTk5qKqqQmZmJtRqNTIyMpCWloa6ujoEBwd32u/NN9/EuXPnkJ6e3vaJ1WqEh4fbr2dnZ+OBBx4AAJjNZpjNZgwcOBABAQGCL4mInO2qUaw2yWh9lEKRCHLSKIxogbraRgtstpZ6MkSkDNl/JuXm5mLs2LFQq1tyntTUVFgsFuTn53fZLzY2FgkJCV22qa2txb59+/Dggw8CaBl9CQoKwqOPPgqNRoMJEybg6tWrXT6G2WxGTU1NmwsR9VytSX5JfaDlQ8WRaRdPt2LaMKH2X5bzM4hISbI/ZcrKyhAd/fWWQpVKhYiICJSXi5Xb7sjHH38MnU5nX+diMBhgMpkwefJknDp1CiqVCkuXLu3yMVatWgWdTme/JCYm9jguIgLSXxZbvJsuuOjVW8wZP0Co/dXqBvCAaSLlCP2ZdOvx1zabzSmHju3cuRMPPPCA/bEGDhyIkpISLFy4EElJSViwYAH27t3b5WMsX74cRqPRfrl8+XKP4yLq7RoaLag2iS3f9cXRF6DlOIRgf7HPu4oax0+pJqKuyf6kiYuLQ2Vlpf26xWKBwWDo8WFMAPDBBx/gW9/6lv26v79/mwW9ERER3U4JBQYGIiwsrM2FiHrmN/8pEGp/10APWfuikE/+Z5JQ++MlRoUiISLZCUx6ejoOHTqE5uaW+fCCggKo1WqMHDmyRwEcPXoUJSUlmDp1qv229evXY9Kkrz8oiouL22yxJiLX+M8RsSni4EDfro0ZGx4k3Ed0DRERySM7gUlLS0NMTAwyMjJw5coVrFy5EjNnzoRWq0VNTY3DJ6ju3LkT48aNQ2Tk1we+TZw4EXl5ecjKysKZM2ewbt06zJs3z6HHJyLHiJbEX3JvskKReJa9SyYKtT/ILdVEipCdwKhUKmzevBlZWVlITk6GyWTCmjVrAAAjRozAjh07HApg586d9t1HrYYNG4Y33ngDixYtQnp6Ou677z5ZNWeIyHnufmWPUPv5kwcrFIlnSerTedmIzlhY2I7I6YTGe0eNGoWjR4+2u72oqKjLfnq9vt0C4Fb79+/v8PZ58+Zx1IXITUSnPQJVAuceyanN4uFnBoUHqWAQWNx8urQGw/vpFIyIqPfx7QlrInLIore/EGp/4Jn7FIrEQ9ySdB1ssuLAucpOGrdXVmNiAkPkZL6535GIeiSnUGzdRkxYoEKReCZHznlqsvCYRyJnYgJDRG1U1TYKtU9LjlAoEs+WPlisYF/+pWqFIiHqnZjAEFEbd//uY6H2b8wZq1Akni1Q7pqfr9wwczs1kTNxDQwR2TU0WtDQJH/HTJBaQkiQAh8jTjqEUWkpfUNwtrxWdnvTAw91P/3k4QuYiTwFR2CIyO6Z/3R9OOutPvv1/QpF4h0SI8S2VOcKLPwloq4xgSEiu21HKoTa67T+CkXiHUSPgrMBaGZNGCKnYAJDRACAyhqzUPsJvXTx7q1Ez386dtmgTCBEvQwTGCICAExb1/WJ77d6rZcu3r2V6PlPVfViu7yIqGNMYIgIAFBeZxFqr8jiXS8VFRwg1J41YYh6jgkMEeFiRZ1Q+3tTohWKxDvdkRAu1D7/UpUygRD1IkxgiAhT1mYLtf/zY6OVCcRLqVViq3lvmMVGu4ioPSYwRIQmwRkNTh+1Fyb4nnAaiahnmMAQ9XKiRwfs/vkEhSLxbiP7i+3K+vS82HlTRNQWExiiXu7+P4gdHTAkPlShSLybv5/Yx2mjxQqrjTVhiBzFBIaoF7NYbbjeIP9LNFrLj4yujEsSqwlz8ZrY4mki+ho/jYh6sf/mXxVq/8EvJikUiW8IFVwHU8QEhshhTGCIerHF7x8Vah8TFqhQJL7D30/+jiQbwGkkIgcxgSHqpWpNzULt79aHKBSJbxmXJFYjh9NIRI5hAkPUSz3xjzyh9n+Zd7dCkfiWIH+xj1UmMESOYQJD1EsdvHRDqD1rv8gnMo0E8IRqIkcwgSHqhYz1TULtv39nP4Ui8U2i00jHeUI1kTAmMES90A/+dkCofeb0OxSKxDeJTiNd5wnVRMKYwBD1QidK64XaawL8FIrEd4lNIgHcjEQkhgkMUS9ztapBqP2r309VJhAfd9dAsWmkcqPYvwtRb8dVeUS9zOQ/7Gl324b3V3Te/kTfluGE7dsVjMr3aAPFRq1OlNYgNlyjUDREvocjMES9jMki1l4SnQshu3jBhIS7kYjkYwJD1IuITh9FBvsrFEnvMKSv2MGXJ64YlAmEyAcxgSHqRab+ca9Q+xEJEQpF0jv4qcSGr67VcTcSkVxMYIh6kdomsSkKteAXMLUXpwsSal9VyySGSA4u4iXqJS5WiJWsvyNe1/aG6dOdGE3vcVtsGEqNJtntv/2nbBz49RQFIyLyDRyBIeolJq3JFmrfJ0xs5IA6JjqNdKVGrEoyUW/FBIaolxCZPFJJ3H3kTAOixHYjNTQKbhUj6oWEEpj8/HykpqZCo9FgypQpqKiokNXv3LlzGDp0KAwGQ5vbY2NjIUmS/TJmzBj7fbt370ZKSgqCg4Mxe/Zs1NeLVQ4loq+Jnn007NbpI+qRQTFiu5Ge3XJMoUiIfIfsBMZqtWLWrFmYNm0aCgsLodFosHjx4m77TZ06FYMHD8bp06fb3WcwGHD69GlUV1ejuroa2dnZ9tsfeeQRLFmyBCdPnkRRURFeeukl+a+KiNqY+7dcofZ9Qjl95EwqweGsfxeUKBQJke+QncDk5OSgqqoKmZmZSEhIQEZGBrZs2YK6uq4XBr755pvYv39/u9vNZjPMZjMGDhyI8PBwhIeHIyQkBACwZcsWJCQkYP78+dDr9Vi2bBk2bdok+NKIqNWxUvkLeP04faSIu5PFjhZobLYqFAmRb5CdwOTm5mLs2LFQq1s2LqWmpsJisSA/P7/LfrGxsUhISGh3u8FgQFBQEB599FFoNBpMmDABV69etT/X+PHj7W3HjRuH4uJiXL58WW64RPQV0d1H+ugQhSLp3TT+YkcL/CX7vEKREPkG2QlMWVkZoqO//gtCpVIhIiIC5eXlDj2xwWCAyWTC5MmTcerUKahUKixdurTD54qKigKALp/LbDajpqamzYWIgMmCu48GRGmVCYSE/OHjs+4OgcijCS3itd1y3rvNZoPk4FjzwIEDUVJSgoULFyIpKQkLFizA3r1fVwm9+blaf+7quVatWgWdTme/JCYmOhQXka8RmYgIDvATXq9B8sUKbk2vNTUrFAmR95OdwMTFxaGystJ+3WKxwGAwIDY21qEn9vf3R1xcnP16RESEfdTk1ue6fv06AHT5XMuXL4fRaLRfON1EJF7VdWxSlEKREADcHhcm1P6pt79QKBIi7yc7gUlPT8ehQ4fQ3NzyF0FBQQHUajVGjhzp0BOvX78ekyZNsl8vLi6GXq+3P1deXp79voMHD0Kv16Nfv36dPl5gYCDCwsLaXIh6u1lvtF9A3xXRomskRvT9zS68rlAkRN5PdgKTlpaGmJgYZGRk4MqVK1i5ciVmzpwJrVaLmpoaNDWJ1ZmYOHEi8vLykJWVhTNnzmDdunWYN28eAGDGjBkoLS3F66+/jqKiIqxevRpz5swRenwiAi5ek1/CnuceuUaIv9j7zGkkoo7JTmBUKhU2b96MrKwsJCcnw2QyYc2aNQCAESNGYMeOHUJPPGzYMLzxxhtYtGgR0tPTcd9999nryuh0Orz33ntYs2YNhg4diqSkJCxfvlzo8Yl6O9FqrkMFpzfIMaOTxLZT//TNzxSKhMi7CR3mOGrUKBw9erTd7UVFRV320+v17RYAA8C8efPsoy63mjJlCgoLC0XCI6KbPPP+YaH2MSxe5xL+fmInuORdNCgTCJGX41lIRD5q27HK7ht9JTokgMXrXCgpKliovcUqcpIVUe/ABIbIB4nuPuofKfaFSj2TFCP2fmd9wV2VRLdiAkPkg2a+mi3UPkIboEwg1CHRWju/+M9xhSIh8l5MYIh8UJFB/q7AID+efeQOdwic+G0Dp5GIbsUEhsjHiO4+mjGy/VllpLw+glV5dx6+olAkRN6JCQyRj8n47zGh9pnfHq5QJNQV0VGvX/xb7N+VyNcxgSHyMVsLSmW31aolaALETkkm5xGpHSg2rkbk+5jAEPmQxmYrGi3y10qs/8FoBaOh7nxTsKid6O4yIl/GBIbIh/x133mh9hOG9FEoEpJDGyg2+vW9N3IVioTI+zCBIfIhr2efk91WBR7e6AnCguQXRD9/rV7BSIi8i9BRAkTkuSxWG+oarbLbp6VEKRgNyTWyfwRyzn5dNXnD+yu6bG85/jr8dmQpHRaRx+MIDJGP+PiY/MW7APDaY2MUioREiJ6N9GVZjUKREHkXJjBEPmLZlvYHrXYlRGDqgpQlMpFXajQpFgeRN2ECQ+QjjGb500eLJw9WMBISNTZJbDpPtFghkS9iAkPkAyprzELtF9ybrFAk5IhQwdGw3wiOthH5IiYwRD5g2rq9stuqVUCAmv/1PU2QwL/JfwSKFRL5Kn6KEfmA8jr5UwozefaRR7pTcBqp1tSsUCRE3oEJDJGXE50+WjGDZx95okDBUbGfvX1YoUiIvAMTGCIv9+1Xc4Ta8+wjzxWukb8WJu/CdQUjIfJ8TGCIvFxpTZPsttw47dmSokNktzULnHlF5IuYwBB5MdF1EOse/oZCkZAzRAYHCrUXnT4k8iVMYIi82M/eEVsH8cCofgpFQs4gSWJF7cb99mPFYiHydExgiLxYztlrsttK4OGN3iBM4y+7rRUsake9F6fEibxUQ6MF3dXevflgwOFxOmD6G8oGRT2Wmhje5nDH7qzcfgKrZnFqkHofjsAQeamV208Kte+rC1IoEnIm0cMdPzldrlAkRJ6NCQyRl9px7KpQe4mzR14jUit/GqmiVv4uNCJfwgSGyAtZrDbUCBzeGKRm7RdvMiIxQqg9dyNRb8QEhsgL7TtTIdT+zqRIhSIhJagFF1t/Z/1+hSIh8lxMYIi80B8+OiPUXrRMPbmfLkj+NNJVI0dgqPfhpxqRFzpVekN223idRsFISCmp/cOF2lfVNioTCJGHYgJD5GVqTc2wClSRHxIbqlwwpBjR3Uiz1u9TKBIiz8Q6MEReZv5bB2W3VYHF63zRzfV92jiw7uuft293TTBEbsIRGCIvc+CCUXbb4CD+jeLNYoIDhNqzKi/1JkIJTH5+PlJTU6HRaDBlyhRUVMjbCXHu3DkMHToUBoOhze2HDx/GqFGjEBYWhhkzZqCqqsp+X2xsLCRJsl/GjBkjEiqRT2pslr91GgAGCZxuTJ5neEK4UPuDF68rEwiRB5KdwFitVsyaNQvTpk1DYWEhNBoNFi9e3G2/qVOnYvDgwTh9+nS7x3vssccwZcoUHDt2DBUVFXj++eft9xsMBpw+fRrV1dWorq5Gdna2/FdF5KP+mnNOqH1UiNjpxuRZ/FQSRGYALSKLo4i8nOwEJicnB1VVVcjMzERCQgIyMjKwZcsW1NXVddnvzTffxP797WsUXLhwAWfPnsVzzz0HvV6PuXPn2tuZzWaYzWYMHDgQ4eHhCA8PR0gI/5IkemPfBdltvzU0mtV3fcBAwVE0JjHUW8hOYHJzczF27Fio1S1z6qmpqbBYLMjPz++yX2xsLBISEtrdrtVqsXbtWgQHBwMArl+/Do2mZbunwWBAUFAQHn30UWg0GkyYMAFXr4qVTSfyNRarDbVm+Wsc1j12p4LRkKv0j9IKtT9TLn+LPZE3k53AlJWVITo6+uuOKhUiIiJQXu7YQWLx8fF4+umnAQBGoxEbNmzA3LlzAbQkMCaTCZMnT8apU6egUqmwdOnSLh/PbDajpqamzYXIl+w9Lb/6rp8EBLB4nU9QSRJEBtIqakyKxULkSYQ+4Ww2W7vrUg/HqOvq6jB9+nSMHDkSCxYsAAAMHDgQJSUlWLhwIZKSkrBgwQLs3bu3y8dZtWoVdDqd/ZKYmNijuIg8zYrtx2W3HR7P2i++JCkmWHbbZk4hUS8hO4GJi4tDZWWl/brFYoHBYEBsbKzDT15fX4+pU6dCq9Xi3XffhUrVEo6/vz/i4uLs7SIiIrodUVm+fDmMRqP9cvnyZYfjIvJElw3yy8U/NKKfgpGQq+mj5CcwANfBUO8gO4FJT0/HoUOH0NzcDAAoKCiAWq3GyJEjHX7yhQsXIjg4GNu2bUNQUJD99vXr12PSpEn268XFxdDr9V0+VmBgIMLCwtpciHyFsb5JqP28u5MUioTcQSU40n2mjFPo5PtkJzBpaWmIiYlBRkYGrly5gpUrV2LmzJnQarWoqalBU5PYB+zRo0exdetWvPrqq2hoaIDBYLDXiZk4cSLy8vKQlZWFM2fOYN26dZg3b57Q4xP5kjl/PSC7bUyIP9e/+CCRisolRq6DId8n+1NOpVJh8+bNyMrKQnJyMkwmE9asWQMAGDFiBHbs2CH0xFu3boXRaERKSgoiIiLsFwAYNmwY3njjDSxatAjp6em47777ZNWcIfJVx8vqZbf9eMm9CkZC7vLNpCih9qKjdkTeRrLdujLXR9TU1ECn08FoNHI6ibxaQ6MFtz+/S3b7opcf+vrK9OkKRETu8vFp+bs+1y/+A7YsSlcwGiJlyP3+5jgzkYdbuf2E7LasW+fbIrT+stueLuE6GPJtTGCIPNzmL67IbjswWqzoGXkXkd1IJp7rSD6OCQyRB2totKBZYJJ38/y7lQuG3C4yWOxsK66DIV+mdncARNS55e93fVTHrSJDAhSKhDyBJLVc5Kxc3PD+ChzZshL3pPTpvNH27c4LjsjFOAJD5MG2HpN/fMAvp6YoGAl5iriwoO4bfaXJYmNRO/JZTGCIfMRP0ge5OwRygSGxYrsqWdSOfBUTGCIPVWaQX4xMAg9v7C38VBLUAkXtynm4I/kofuIReajJa7o+wPRmc+/i4aW9SaxO/jSShTNI5KOYwBB5IIvVhrpGq+z2U4fGKxgNeZrBfcROG29o4p5q8j1MYIg80N7T8hfvAsA3B4mVmSfvJnIuEgB8duGaQpEQuQ8TGCIP9MKOk7Lbhgf5CX+hkfcT+Sdvlj+YR+Q1mMAQeaBSQ4Pstj+9J1nBSMhTfTMpWqi9mVkM+RgmMEQexmK1QWD5C55IH6hcMOSxtIF+Qu0PXbyuUCRE7sEEhsjDfFhQIrttcIAft0/3Yhp/+UmMiSMw5GP4yUfkYRZsPiK77V0DI5ULhDxeSl+x3UjNrMpLPoQJDJEHET18b+3sUQpFQt4gOkTscMcTVwzKBELkBkxgiDzInP/NE2ofEsTzWHszSQL8JPnbkarqGxWMhsi1mMAQeZDjV2tlt03po1UwEvIWfQWq8nIGiXwJExgiDyFy9hEALJ86VKFIyJsMEVwHw+3U5CuYwBB5iAf/tE+o/YTb+ygUCXkT0SKGBy9wOzX5BiYwRB6iSmABb6BK/IuLfFeQwFb6JosVNk4lkQ9gAkPkASyCixPuTharwkq+7c4ksbOwqurMCkVC5DpMYIg8wM78K0Lt//TYaIUiIW8UKFjM8OK1OoUiIXIdJjBEHuDp948Jtef2abqVTuMvu62hQazeEJEnYgJD5GYWqw0i+0LuiA9RLBbyXqmJ4ULtuQ6GvB0TGCI32/dlhVD7TU+MVygS8mb+fmIf51wHQ96OCQyRm/3+ozNC7XVa+VMF1LvEhckvaneipEbBSIiUx4l0Ijc7WXpDdtsl9yZ/fWX6dAWiIW92W1wYSmvkFURssljR0GiBJkD+idZEnoQJDJEbdXZ444b3V3R4+6QTfYBXWf+FOuankiABkLu85dktR/GH7/NAUPJOnEIicqMfbpB/eKMEQCVwcB/1TmqBAof/LihVMBIiZTGBIXKjoyXyD2/UcqifZBBdI9XZKCCRp2MCQ+QmDY0WofajB0QqFAn5kuH9woXaz/v7QWUCIVIYExgiN8ncdlx2WwlAgGC1Veqd1F+tg5Gr4Ap3I5F3EvpEzM/PR2pqKjQaDaZMmYKKCnn1K86dO4ehQ4fCYDC0uX337t1ISUlBcHAwZs+ejfr6eln3EfmCXSfkrz8YkRCuXCDkc+7opxNq39gsUkqRyDPITmCsVitmzZqFadOmobCwEBqNBosXL+6239SpUzF48GCcPn26ze0GgwGPPPIIlixZgpMnT6KoqAgvvfRSt/cR+QqjWf6XRnRIoIKRkK+JCZVfDwYA/pJ9XqFIiJQjO4HJyclBVVUVMjMzkZCQgIyMDGzZsgV1dV0fCvbmm29i//797W7fsmULEhISMH/+fOj1eixbtgybNm3q9j4iX3CxQv5hehIAbj4iEaK/L2s/OatMIEQKkp3A5ObmYuzYsVCrW0rHpKamwmKxID8/v8t+sbGxSEhI6PDxxo//uiT6uHHjUFxcjMuXL3d5H5EvmLwmW3Zb0ZOGiURZbOKLyoncTfYnY1lZGaKjo7/uqFIhIiIC5eXlDj3xrY8XFRUFACgvL+/yvs6YzWbU1NS0uRB5ItHDG/tFaBWLhXzXuKQoofYrtp9UKBIiZQhV4rXdcnypzWaD1IOx7Zsfr/Xn1sfr6r6OrFq1CitWdFy9lMiTVN07BRuuGmW3HxDFBIbEhQaJFVrfdbwEL88aoVA0RM4newQmLi4OlZWV9usWiwUGgwGxsbEOPfGtj3f9+nUALVNOXd3XmeXLl8NoNNovnG4iT3WqTGx0kNV3yVHRIQGy2xpMnEIi7yI7gUlPT8ehQ4fQ3NwMACgoKIBarcbIkSMdeuL09HTk5X1dRv3gwYPQ6/Xo169fl/d1JjAwEGFhYW0uRJ7GYrWhySL3pBogTPCvaKKbiRa1qzU1KxMIkQJkJzBpaWmIiYlBRkYGrly5gpUrV2LmzJnQarWoqalBU5NYOeoZM2agtLQUr7/+OoqKirB69WrMmTOn2/uIvNl/DouNDI7sH6FQJNQbiJyLBAA/feuQQpEQOZ/sBEalUmHz5s3IyspCcnIyTCYT1qxZAwAYMWIEduzYIfTEOp0O7733HtasWYOhQ4ciKSkJy5cv7/Y+Im+27N/yq+8CgL8fdyBRz0QFy59GyrtQpWAkRM4lND49atQoHD16tN3tRUVFXfbT6/XtFgADwJQpU1BYWNhhn67uI/JW8iePWLyOnOOOhHBkn5FXNR1oOdxR9EBIInfgn3dELiK6vmC4YDl4oo6ITiP96B883JG8AxMYIhf5+TuHhdqLfvEQdSY00E922/zLrKFF3oEJDJGL7Dl7TXbb4AD5XzhE3Rk1IFKovbFebFMGkTswgSFyAdHpo9GCXzhEXRFdDD6P00jkBZjAELnAU29/IdQ+gOcfkZP5C0xJnrjKaSTyfPyUJHKBnMLrstty+oiUMHag/LORmqwtRReJPBkTGCKFNTZbhbZPc/qIlKDxF0uMdx6+olAkRM7BOuVECvv7gQtC7Tl9REqJDQtCWY3Jfn3D+10cgPs+gNv7Atu3Kx8YkQP4SUmksPV7zspuy63TpKQ4nUaovbnZqlAkRD3HBIZIQY3NVtxolD+BNE5gnQKRqEiBYwUA4NBF+Wu3iFyNCQyRgv5XcPpIdJ0CkQhJAoIEpihNHIEhD8YEhkhBf9h9RnZbnn1ErhAdKvZ71tBoUSgSop5hAkOkkFpTM5oFth/x7CNyhcF9QoXaP7f1mEKREPUMExgihSx8R6x4HRfwkiv4qSShWkNbC0oUjIbIcUxgiBSy76z8BZBarn0hFxJZLN5sY1E78kxMYIgUIHr2keiwPlFPqCSx0b6DF7gbiTwPExgiBfzi3Xyh9qILK4lcadUHp9wdAlE7TGCIFPDZxSrZbdUqCYJ/EBP1WN8w+Unzias30Mgt1eRhmMAQOZnFasMNs/ytp33CghSMhqhjQ+PEdr39fb9YTSMipTGBIXKyvHPXhNoP6cv1L+R6foK73l7LOadQJESOYQJD5GRL3z8qu61aJQl/kRA5i8h26hqThdNI5FGYwBA5UUOjBeU1ZtntWbyO3Gn0gEih9v974LxCkRCJYwJD5ESiVUujgrn7iNwnQOBcJABY+5H8k9WJlMYEhsiJth2RX7U0XBPA3Ufkdt8QGAU0W3g2EnkOJjBETmKx2tAksEQgKSZYuWCIZIoOFdsFt3L7SYUiIRLDBIbISQ6cqRRqH6kNUCgSIvkkCQj0k/9V8MnpCgWjIZKPCQyRk/zkn4dktw0NVHP6iDzGbXFhsttW1MpfpE6kJCYwRE5grG9Co8D0UayOxevIc0SHiC0m5zoY8gRMYIic4Ef/+EyofWKkVqFIiMRJEuAvMI1098ufKBgNkTxqdwdA5AtOXDXKbhuh9Rc+DZhIad8cGIX9hfLWcVXVN6HW1IyQoK++QqZP777T9u09iI6oPY7AEPVQY7NVaPpoZP8I5YIhclCgYE2YhZsOKxQJkTxMYIh66Nf/OS67bWK4hqMv5LFEkph9gmd+ETkbExiiHrBYbXg//4rs9j8cr1cuGKIeEl3MW2tqVigSou4JJTD5+flITU2FRqPBlClTUFHRfT2A3bt3IyUlBcHBwZg9ezbq6+sBANnZ2ZAkqc1Fr9fb+8XGxra5b8yYMWKvjMgFcgRrv/w/JjDkwVIET0Z/+t18hSIh6p7sBMZqtWLWrFmYNm0aCgsLodFosHjx4i77GAwGPPLII1iyZAlOnjyJoqIivPTSSwCAtLQ0VFdX2y8vvvgi4uLi2vQ9ffq0/f7s7GzHXiGRgl7ZJb8qqUoSP3uGyJX8VBLUAqej7/1SLIEncibZn6Y5OTmoqqpCZmYmEhISkJGRgS1btqCurq7TPlu2bEFCQgLmz58PvV6PZcuWYdOmTQAAtVqN8PBw+yU7OxsPPPAAAMBsNsNsNmPgwIH2+0NCQnr4Uomc73xFvey2kVp/BSMhco67k6Nlt7WC00jkPrITmNzcXIwdOxZqdcu2udTUVFgsFuTndz6EmJubi/Hjx9uvjxs3DsXFxbh8+XKbdrW1tdi3bx8efPBBAC2jL0FBQXj00Ueh0WgwYcIEXL16VeiFEblCk01+26yfTVAuECIn8fdTQWSZ+ew3chWLhagrshOYsrIyREd/nZmrVCpERESgvLxcdp+oqCgAaNfn448/hk6ns69zMRgMMJlMmDx5Mk6dOgWVSoWlS5d2GZ/ZbEZNTU2bC5GSth+Vf/I0AMSGs/oueYd+4RrZbU+U1cImkMgTOYvQhLztlt9Sm80GqZstoTf3af351j47d+7EAw88YL994MCBKCkpwcKFC5GUlIQFCxZg7969XT7PqlWroNPp7JfExETZr4tIlMVqw9PvFshu/79zRysYDZFzDRZczFtVx/ORyPVkJzBxcXGorPx6wZbFYoHBYEBsbKzsPtevXweAdn0++OADfOtb37Jf9/f3b7OgNyIiotsRleXLl8NoNNovt05TETlTzplKWAX+6px4e1/lgiFyMj+VJDSNVFhRq1gsRJ2RncCkp6fj0KFDaG5uWbBVUFAAtVqNkSNHdtknLy/Pfv3gwYPQ6/Xo16+f/bajR4+ipKQEU6dOtd+2fv16TJo0yX69uLi4zRbrjgQGBiIsLKzNhUgpz205JrvtWH0k/AR2dhB5Am2An+y2teZmTiORy8lOYNLS0hATE4OMjAxcuXIFK1euxMyZM6HValFTU4OmpqZ2fWbMmIHS0lK8/vrrKCoqwurVqzFnzpw2bXbu3Ilx48YhMjLSftvEiRORl5eHrKwsnDlzBuvWrcO8efMcf5VETmSx2nC1Rv6Q+c8mJSsYDZEyRg+I7L7RTSpqTApFQtQx2QmMSqXC5s2bkZWVheTkZJhMJqxZswYAMGLECOzYsaNdH51Oh/feew9r1qzB0KFDkZSUhOXLl7dps3PnTvvuo1bDhg3DG2+8gUWLFiE9PR333XdftzVniFwlr1CshPp4gW2pRJ4iQK0SOvbiRIn8A02JnEGy3boy10fU1NRAp9PBaDRyOomcavhzO1DbfsCxQwnhQTjwzOS2N8o5uZfIA1y6Xo/Cihuy20+6rU/nSQ9PoyaZ5H5/sywokYBaU7Ps5AUAXvj2cOWCIVJYYqT87dQAcLlKfmFHop5iAkMk4IG12ULtJ9zWR5lAiFxAJUnwE1h/XswEhlyICQyRTI3NVlwxyF+8G6SWuPuIvF5StPxjXMzNVlh9c1UCeSAmMEQy/SP3glD7T5bcq1AkRK7TP0or1L7oWufn4xE5k9rdARB5i6FPzsGGZmuXbZ54OMP+cz/B9QNEnkglSYjQ+qO6Xt7irwvX6jAwhofvkvI4AkMkQ2OzFeZukpebpSaIlWIn8mQj+0cItW9osigUCdHXmMAQyfBm3kWx9o/fpVAkRK6nksSOFvj0vFitJCJHMIEhkmH1ri9ltw3wA3RafwWjIXK94f10sttabeBiXlIcExiibhjrm9Aof/YIf/vhncoFQ+QmfUKDhNoXXediXlIWExiibqS9/LFY+8ExCkVC5D6SBIQGyj/g8UIlExhSFhMYoi40NFpwQ2D4JSlKw9ov5LNidWJbqpssAkOXRIKYwBB14dktx4Tab12UrlAkRO4nerTA/sJKhSIhYh0Yoi79u6BEqL3u+zMVioTI/VSShAA/FRpljqxYbUCz1Qa1SpJ3iCkPfCQBHIEh6oRRZuGuVlHceUS9wADByrwnrhqUCYR6PY7AEAEd/nV4tqgKGxrkJzF3JIoV+yLyRomRWhRW1Mpuf622UcFoqDfjCAxRJ2oEkhc/CS3D5EQ+TiVJiAkJFOojUsWaSC4mMEQdsFhtECnDlZ7SR7FYiDzNiIRwofaHLlYpEwj1akxgiDpw7IpBdlsVR1+ol5EkICxI/goEUzPPRiLnYwJDdAubDbheJ3/ePj6cp05T7zNqQKRQe04jkbMxgSG6RXW92KLDwX148jT1PmqVJFS0Mfcca8KQczGBIbqFyPRRoFrFyrvUa40QPODRYuUBj+Q8TGCIbtJksaJZ4EN2aLz8D3AiXxMZLLYb6Uz5DYUiod6ICQzRTURLn0dqAxSKhMjzSRIgMv5YYmhQLBbqfZjAEH2lyWKFyAh3WJA/JM4eUS83flC0UHse8EjOwgSG6CtHLhuE2o8awMq7RJoAP6H2By9cVygS6m2YwBB9xShQeVcCa78QtYoOlj+Vam62wmrjYl7qOSYwRBAf1hbZfUHk64YLVubNPXdNmUCoV2ECQwQgV3DxbnRokEKREHkftUqCv8CIpLlZbLcfUUeYwFCv19BoQbPAZ2mcLoiLd4luIToKc+KqQZE4qPdgAkO93rfXZQu1vz0uTJlAiLyYaEmBa7ViFa+JbsUEhnq1xmYrCq+bZLdXSYCKwy9E7UgS0DdMbGrV1MQt1eQ4JjDUq/157xmh9uMHxSgUCZH3GxYvNjr56Xku5iXHMYGhXu3VTy4ItQ/y538Zos6oJAl+AgOUFpuNW6rJYUKfxvn5+UhNTYVGo8GUKVNQUVHRbZ/du3cjJSUFwcHBmD17Nurr6+33xcbGQpIk+2XMmDGy+hE5Q2OzFSIfnd/g1mmibo1PFhulvHSdn+3kGNkJjNVqxaxZszBt2jQUFhZCo9Fg8eLFXfYxGAx45JFHsGTJEpw8eRJFRUV46aWX2tx/+vRpVFdXo7q6GtnZ2bL6ETnDhN/tEWrPrdNE3QtUi41Snq+sVSgS8nWyf9NycnJQVVWFzMxMJCQkICMjA1u2bEFdXV2nfbZs2YKEhATMnz8fer0ey5Ytw6ZNmwAAZrMZZrMZAwcORHh4OMLDwxESEtJtPyJnqDU1o+yGWXZ7jb+KW6eJZBqZGC7UPo9rYcgBshOY3NxcjB07Fmq1GgCQmpoKi8WC/Pz8LvuMHz/efn3cuHEoLi7G5cuXYTAYEBQUhEcffRQajQYTJkzA1atXu+1H5AwTfvexUPuxSVEKRULkeyKDA4Xa1zdaeMgjCZOdwJSVlSE6+utTR1UqFSIiIlBeXi67T1RUy5dAeXk5DAYDTCYTJk+ejFOnTkGlUmHp0qXd9uuM2WxGTU1NmwtRRxoaLahqsMhur1ZJ8Pfj4l0iuSQJGBClFepzqKhKoWjIVwl9KttuWS1us9kgdTOufnOf1p8lScLAgQNRUlKChQsXIikpCQsWLMDevXu77deZVatWQafT2S+JiYnyXxj1KjP+nCPUPn0wt04TiUqOCRVqX99ogYXHC5AA2QlMXFwcKiu/Pi/GYrHAYDAgNjZWdp/r11uOUY+NjYW/vz/i4uLs90VERNhHTbrq15nly5fDaDTaL5xuoo40NltxtrJBdns/SYIfT50mEiZJQJDaT6hPXiHXwpB8shOY9PR0HDp0CM3NzQCAgoICqNVqjBw5sss+eXl59usHDx6EXq9Hv379sH79ekyaNMl+X3FxMfR6fbf9OhMYGIiwsLA2F6JbvZlXJNT+7uTo7hsRUYfGJkUKtX/uv8cUioR8kewEJi0tDTExMcjIyMCVK1ewcuVKzJw5E1qtFjU1NWhqamrXZ8aMGSgtLcXrr7+OoqIirF69GnPmzAEATJw4EXl5ecjKysKZM2ewbt06zJs3r9t+RD2x/ehVofYBgltCiehrAWqV0P+housmNDZzMS/JI/s3S6VSYfPmzcjKykJycjJMJhPWrFkDABgxYgR27NjRro9Op8N7772HNWvWYOjQoUhKSsLy5csBAMOGDcMbb7yBRYsWIT09Hffdd5+9rkxX/YgcZbHacPyq/MXd0cFih9MRUXsTBNeQpa7YrVAk5Gsk260rc31ETU0NdDodjEYjp5MIADDztVzkFxs6vG/D+yva3TZxSB+ouf6FqMcOFVXB2NB+lP5WTzycAQDIf/Z+RIbwD4jeSu73t9qFMRG5x/TpsFhtWHim+6MvWgUH+DF5IXKS1MRw5Jyt7L7hV8b99iMU/vYhBSMiX8AJfuoVcs+J7W745kAu3iVyFn8/ldAfBE3WlmrZRF1hAkM+r9lqQ6NAlc+B0SE8NoDIyUTrKT382n6FIiFfwQSGfF7eOflD1wCQFB2sUCREvZdoPaUvK+pZ2I66xDUw5NMaGi1otMj/EIwMDuDoC5FCkiK1uFhVL7v9gbOVuOe2Pm1vnD69+47btwtGRt6IIzDk0yat/kSo/TcSwpUJhIiQ1CdEqP3P/9X5YcFETGDIZzU0WlB6o/utm638JPFhbiKSTyVJ8PeT/3/M2GBB1pESBSMib8YEhnzWE29+JtT+7mQe2kiktPGDxHb4PfVuAdfCUIeYwJBPslhtyD1fLbu9SuKxAUSu4O+nguhA5+Q/7FUmGPJq/MQmn5T8651C7Sek9Om+ERE5heiW6qLrDawLQ+1wFxJ5Lgd3G8z9636IDDiLFtkiop7x91NB4++HhiaL7D6PvL4fOxffq2BU5G04AkM+paHRgv0X5B/YCABpyay6S+RqdydHC00lnSpnXRhqiwkM+ZQ7MncJtddp/LnziMhNvpEYIdT+T5+cVSgS8kZMYMhnVNaY0Sz/xAAAwJgBkcoEQ0TditSKnTi97pNzsHEQhr7CNTDkM2asFzs7JSTAj1V3idxIkoBIrT+q6lvqNW14f0W3fbJVEu4dwkX3xBEY8iElRrNQ+zFJUQpFQkRyiU4jWaw2mEWHWsknMYEhn5B15KpQe7VK4s4jIg/gp5IQExIo1OfghWsKRUPehAkMeT2L1Yan3j0i1CdNsA4FESnnG4nhCBQoJNlksXEtDHENDHm56dORX1SFDQ3yzzwKUvtx9IXIwwyL1yG/WH717L1nyjHptr4KRkSejiMw5NUsVhuMAskLAIxP5toXIk8ToQ0QqgtjtQGny4zKBUQejwkMebXcc2Jz4RFaf6i49YjI40gSMDReJ9TnarUJVs4l9VpMYMhrnS2/gUaL2G6Ekf3FdjwQkevEhgUhLMhfqM+eLysUioY8HRMY8kpWmw3FVfVCfRIjtRx9IfJwY5PEi0t+duG6ApGQp2MCQ16poNgg1N7fT4UhfUOVCYaInCpepxFqf8PcjGaek9TrcBcSuYeck6Y7UV5jQnV9o1Cf9ME8sJHIWwyJDUWJsUGoT965a5iQwvIIvQlHYMir2GzA8atiOw/6hgZy6ojIi/ipJESHiJ2T1Gix8rTqXoYJDHmVY1cMwn2G9wt3ehxEpKzUxAjhPzz2neWC3t6ECQx5DavNhspasfOOhsfreGAjkZe6R3BKyGIDSo0mhaIhT8MEhrxGzplKofZBaj/E6oIUioaIlOankhAdLDaVdLLEyKmkXoIJDHmFUkMDLIIFq1hxl8j7pTpQu2nkil0KREKehgkMeTybDThZWiPUJyFCw4W7RD5iguDhqzVmKya88olC0ZCnYAJDHu+TL8uF+9wWG6ZAJETkDgFqFQL8xL6uiqtMMNaLnZNG3oUJDHm0gw5U2ExLZi0IIl/jSI2XMS9+qEAk5CmEEpj8/HykpqZCo9FgypQpqKjofsva7t27kZKSguDgYMyePRv19V+Xfz98+DBGjRqFsLAwzJgxA1VVVfb7YmNjIUmS/TJmzBiRUMmdpk/v/iJDs9WGWnOz0FOrAAT5My8n8kUTh/QRat9kBX688ZBC0ZC7yf6kt1qtmDVrFqZNm4bCwkJoNBosXry4yz4GgwGPPPIIlixZgpMnT6KoqAgvvfSS/fEee+wxTJkyBceOHUNFRQWef/75Nn1Pnz6N6upqVFdXIzs727FXSF4r54x4TYdJt/dVIBIi8gRqlQQ/wbVtn3xZgYZGi0IRkTvJTmBycnJQVVWFzMxMJCQkICMjA1u2bEFdXV2nfbZs2YKEhATMnz8fer0ey5Ytw6ZNmwAAFy5cwNmzZ/Hcc89Br9dj7ty52L9/PwDAbDbDbDZj4MCBCA8PR3h4OEJCQnr4Usmb5J6rhOhGyHsF/zojIu+T5sCxIHdkcleSL5KdwOTm5mLs2LFQq1uOT0pNTYXFYkF+fn6XfcaPH2+/Pm7cOBQXF+Py5cvQarVYu3YtgoODAQDXr1+HRtNygJfBYEBQUBAeffRRaDQaTJgwAVevXnXoBZL3KTWa0NBkFeoTGRwAPxV3HRH5On8/FTT+fkJ9mq3A8Oc/UCgichfZCUxZWRmio7/OfFUqFSIiIlBe3vkOkVv7REW11OUoLy9HfHw8nn76aQCA0WjEhg0bMHfuXAAtCYzJZMLkyZNx6tQpqFQqLF26tMv4zGYzampq2lzI+9hsLYWoRI1yoFYEEXmnu5OjIfrnSm2jFQ+ty1YiHHITodWOtlsKidlsNkjdzEfe3Kf155v71NXVYfr06Rg5ciQWLFgAABg4cCBKSkqwcOFCJCUlYcGCBdi7d2+Xz7Nq1SrodDr7JTExUeSlkYcorLgh3Ed0YR8Reb97bxP/f3+ytA61JrGNAeS5ZCcwcXFxqKz8upS7xWKBwWBAbGys7D7Xr7dsiW3tU19fj6lTp0Kr1eLdd9+FStUSjr+/P+Li4uz9IiIiuh1RWb58OYxGo/1y+fJluS+NPMTnF6tQXFXffcObhAapoebUEVGvo5Ik9I/UCve7a9XHCkRD7iA7gUlPT8ehQ4fQ3NySvRYUFECtVmPkyJFd9snLy7NfP3jwIPR6Pfr16wcAWLhwIYKDg7Ft2zYEBX19Zs369esxadIk+/Xi4mLo9fou4wsMDERYWFibC3mP3PPXUGMSKzqlkoBxSTwugKi3SukbKrwe5obZgq1f8A9cXyA7gUlLS0NMTAwyMjJw5coVrFy5EjNnzoRWq0VNTQ2amtp/+cyYMQOlpaV4/fXXUVRUhNWrV2POnDkAgKNHj2Lr1q149dVX0dDQAIPBAIPBAACYOHEi8vLykJWVhTNnzmDdunWYN2+eU14weZ6z5Tcc2uZ4Twqnjoh6u7uTo4Ursv7i/WPYdaJUkXjIdSTbrQtbupCfn48f/ehHOHPmDCZMmIC3334bMTEx0Ov1WLt2Lb7zne+06/Phhx9i0aJFuHr1Kr797W/j73//O7RaLVasWIHMzMx27VvD2bhxIzIyMtDQ0IAf/OAHeOWVV+Dv7y/7hdXU1ECn08FoNHI0xtVkFqoDAKvNhj1fitd7iQoOwEgu3CUiOPY58tPvZaDwpQe5e9EDyf3+FkpgvAkTGDcSSGA+OV0uXO9FJUmY5MACPiLyXYXlN3BJYA3dEw9nIDRAheMrv6VgVOQIud/frLlObnOgULxYnZ+KyQsRtTe4byiCAwTXwzRacddLLHLnrZjAkFtkn6mAqVmsWB0A3OPAgW5E1DuMGyi+qL/0hgVpqz5RIBpSGhMYcrnPLlxDs1V85jIxUguV4DkoRNR7qCQJA6LEt1ZfMZow6gWeXO1tmMCQSx0prsYNs/iOowC1CkP6hioQERH5ksF9Qh2qD1NV14T0lzkS402YwJDL5J67hmt1jcL9gvxVmDCYU0dEJE9K31BEhwQI97tsMGHePw4qEBEpgQkMuUTO2Qo0NImPvEgA0pKZvBCRmNTECOEidwCQfeY6Mv57QoGIyNmYwJDi9nxZjiaLY7v1HTnvhIgIaCly5+9AnZc38y7hxxs/UyAicia1uwMgDyKnfsv27UIPmXO2Ag6s1wUA9OeiXSLqoQkpffDJl+Xtbt/w/opu+86zvYqNPxqrRFjkBByBIcV8dvG6wyMvoYFqpHDRLhH1kCQBIxJ0DvXNPlOJaev2OTkichYmMKSIfWcrcMPBY+v9JMmheg5ERB3pExqEO/o5lsScKL2BO57/ABZHh5JJMUxgyOlyz11Do4MjL36SxHUvROR0fcMcT2JuNFqR/OudPADSwzCBIad6YG22Q7uNACBS68/khYgU0zcsyKFCdwBgAzB/Uz7+m3/FuUGRw5jAkFNYrDYM+c0OfFlW51D/AD8VRg2IdHJURERtDe4Tijv66RzeIPDz947ixxsPOTkqcgR3IZGYDnYqldeYcPyqEesdfMgAPxUm8IwjInKRvmFB6BMahNzz12ByYMT4ky8rMO1P+5H183QFoiO5OAJDDrPabDh8qRrHrxodfozo4AAmL0TkcpIEjB/k+GaBEyU1SM3chYZGx6bMqeeYwJBDCituYM+XFaiuFz8aoFVipBap/SOcGBURkXyOHv7YymCy4Pbnd+Hxf3zuxKhILiYwJMRiteGzC9dx6Xp9jx4nOiSAhzMSkdsN7hPaoyQGAPacqcTolbudFBHJxQSGZDty2YC9Zypww+xYfZdW2gA/pCZy5IWIPMPgPqGYdFsf+Ps5Xvn7en0zhjy7k1NKLiTZbDafrM5TU1MDnU4Ho9GIsLAwd4fjHTo5SsBqsyHnbKVTCjn1j9Sywi4ReaxSowmnSoyQ+2n3xMMZ7W5L6RuMrJ9NQICaYwSOkPv9zV1I1KnGZis+PX8NTU5IXFQAJt7Wh2cbEZFHi9MFoW9YILLPOH6O29nyOqQ8+wFujw3BfxamQRMgfio2dY8jMPS1m0Zgcs5WoslidcrD+qkk3DuEBeqIyHtU3DDh2BXHd1jeLPOJVfhoyUQmMjLJ/f7m+Ba1YbMBH58ud1ryEh0cwOSFiLxOn9AgjEjQwRljxlcMJtz+/C488SZ3KzkTExiCxWrD/rOVOHypusNj5x0R+VXiwm3SROSt+oQGYdJtfRGu8XfK4318uhJ3rfoYjc3O+QOxt+MamF5u14lS/M97R1HXaMGGHtR0udmASC0Gc6EuEfkASQLG6CNhsdqQc9bxdTGtSo1mpDz7AaaNiMO62SPhp+K6QEcxgellLFYbDl64jk/PX8f5yhv44IRzRlwAQALQP0qLwX2YvBCRb/FTSZh0W1/sO1uJRidMsWcdK8WHJ8uw6N5kPDVpMBMZBzCB6QUam6343wMX8M+8SyitMcneHijC30+F9MHR3GVERD5tQkoMjlw24FqtuceP1Wix4Y8fF+Jv+y8gfXAMBsUE466B0fjmoCgmNDIwgfFRFqsNB89fx+8//BIFl52zkr4zYUFqjE1y/EwRIiJvkpoYbh/NbnDgMMhb1Zot+OBEGQDg1b3n4a+SEBcehP4RWvwkfSDSUmKY0HSACYwP2nmsBMv+fQy1ZmUrQmr8/TBuYBTU/I9FRL2Mn0rC3cnRTk1kWjVZbSiuakBxVQMOnL8OPxXws3sH42eTOdV0MyYwXspiteHzi1WouGFCpCYAX5bfwOXqehSW1+DTC9Xt2m94f4XTnjs4oCVx4XQREfV2rYlMeY0JJ0qMUKKymsUKrP2kZarpdzNHoLTGhENFVdAG+GHWyASMHxzdKxMbJjBeaNeJUqzYfgqlRpPLn/veIX165X8UIqKu9A0LQp/QIFTVmXH8qhHNgtuV5PyR+cTDGXjq3YI2t209UgJtgB9+//AIRAQHouKGCX1CgzA2KdLnP6uZwHiQm0dVokMCARtwrc7c5pdx14lSLNiUr8hC3M6oVRK+OTAaQf4sG0RE1BlJAqJCAjFxSB+U1ZhwqqQGVhcUu69vtGDhO20Tm3CNP+beNQBhQf64XF2PAZFazL1L71PnMzGBcYHGZiv++WkRLlV1/Etksdrw6p5z+EfuRRgamjp8jDhdEJ57aChe2HHKZclL37AgDI/XgTNFRERiYsOC0Dc0CCdKjE6dwpfL0NCEP+851+a2F3ecxrQRsbhvaKxPjNLwLCQBN4+QdPuPP306LFYbvrhUhRum5jZ32eulfJaNXSdK8cx/jsNQ33HicnOfzv6hnPWfI0IbALWfhHCNPxIjtVzjQkTkBFabDV+W3kBZjalHIzIdnXx9K7lTUUDLH8YZ04figeFxbe6/9btu9IAIHL5U3el1ZydCPI3ayTpad9LZPz4AHL1iQOWNjusE2ABcul6PV946hI9OVch6fiWzTH8/FW6PC0Wf0CAFn4WIqHdSSRKGxofh9rgwXLxeh8tV9U47b64nyowmLNiUj9fnjLJ/j3X0XaeS0KYC8a3Xu/ouVJLQCEx+fj4ef/xxnDlzBunp6di0aRP69On6oL7du3fjZz/7Ga5evYrp06fj73//O7RabY/uk8OZIzCdrTtpzTdv/scHgJ+8dQjfX7moR8/ZylkZd0e0AX64LTYMEdoAThMREbmIzQZU1zfiWq0ZJYYG4QW/PXXz94oEIFYXhAO/moSPTpXZv+tERnI6+y50lNNPo7ZarZg1axamTZuGwsJCaDQaLF68uMs+BoMBjzzyCJYsWYKTJ0+iqKgIL730Uo/uczWL1YYV2zted9J624rtp2D56hewodEie1TF1SQAgWoVEsI1uHdIH4wfFI3IYCYvRESuJEktB96m9A3FPSl9MKp/BJKig6EN8HN5LDYApUYTDl643ul3nZzHANp+F7qC7AQmJycHVVVVyMzMREJCAjIyMrBlyxbU1dV12mfLli1ISEjA/PnzodfrsWzZMmzatKlH97na5xerutyu3PqP//nFKgDAb3eecklcIjmHWiVhYHQIJt3WF+mDY3BbXJhXL9wiIvIVrcnMoJgQjB8UjXuH9EFUcIDQZ7wzfHr+eo9Kc9z6XegKstfA5ObmYuzYsVCrW7qkpqbCYrEgPz8f6enpnfYZP368/fq4ceNQXFyMy5cvO3xfYmJih89lNpthNn+95sRobCmfX1NTI/cldqiotBJWc32b2/687eV27aL3h6EmNBBTrhgwpr4Jnad1Yta997z959YRlEExIThfWQvTV0ey3/pcEoCY0EBEhwTC30+FCI0/JAm40dz1QmEiInK/QXEhGGgDqhuaYKhvBACEawNQXWfG5eqGHq+JvPU7DQBMdTfa3F5n7b6ycEePU1RaiWEx/j2Kr/V7u7sVLrITmLKyMkRHR9uvq1QqREREoLy889OMy8rKcMcdd9ivR0W1nJdTXl7u8H2dJTCrVq3CihXt5+w6a98T3+noxotOf5rOFbrwuYiIyLesfaTdTc+tbXv9Ow4+zmNr2zdz1I0bN6DT6Tq9X2gX0q3ZkM1mg9TNAoqb+7T+3NrH0fs6snz5cixZssR+3Wq1oqqqClFRUd3G2J2amhokJibi8uXLTtuS7Y34PrTg+8D3oBXfhxZ8H/getHLG+2Cz2XDjxg3Ex8d32U52AhMXF4fTp0/br1ssFhgMBsTGxnbZp7Ky0n79+vXrAIDY2FiH7+tMYGAgAgMD29wWHh4u45XJFxYW1qt/MVvxfWjB94HvQSu+Dy34PvA9aNXT96GrkZdWshfxpqen49ChQ2hubinKVlBQALVajZEjR3bZJy8vz3794MGD0Ov16Nevn8P3EREREclOYNLS0hATE4OMjAxcuXIFK1euxMyZM6HValFTU4OmpvYLRGfMmIHS0lK8/vrrKCoqwurVqzFnzpwe3UdEREQEm4DDhw/bRowYYQsMDLTdf//9toqKCpvNZrMNGDDAtmXLlg777N6925acnGzTaDS273//+7a6uroe3+dqJpPJlpGRYTOZTG6LwRPwfWjB94HvQSu+Dy34PvA9aOXK98Fnz0IiIiIi3+U752oTERFRr8EEhoiIiLwOExgiIiLyOkxgnOCtt96CJEkoKipydygud/jwYYwaNQphYWGYMWMGqqpcdw6Gq+Xn5yM1NRUajQZTpkxBRYVnHtqptAsXLuCee+5BaGgoJk6ciEuXLrk7JLfat28fJElCdna2u0Nxi6amJjz55JMIDQ3F0KFD8fnnn7s7JJc7fvw47rrrLoSGhmLq1KkoLi52d0guU1painvuuQdHjhyx3+aqz0omMD1kMBjwy1/+0t1huIXVasVjjz2GKVOm4NixY6ioqMDzzz/ffUcv5Mhp7L7qpz/9Kfr3748TJ04gKioKixYtcndIbtPU1ISFCxe6Owy3+v3vf4+ioiIUFBRg9uzZvbLkxXe/+11MmzYNZ86cgV6vx+OPP+7ukFziySefRHx8PPbt22e/zaWflYrvc/JxCxYssH3nO9+xAbBdvHjR3eG4VGFhoQ2Arba21maz2Wzr16+3jRgxws1RKWPPnj22sLAwW1NTk81maykpoNFo7K+9tzCbzTZJkmwnT5602Ww2244dO2xhYWFujsp9fve739nuvfdem06ns+3du9fd4bjFoEGDbEeOHLHZbDbbjRs3bJs3b7ZZLBY3R+U6FRUVNgC20tJSm81ms+Xl5dm0Wq2bo3KNyspK28WLF20AbAUFBTabzbWflRyB6YHDhw9j8+bNeOWVV9wdiltotVqsXbsWwcHBAFqOfNBoNG6OShldncbemzQ1NeGVV15BUlISAN/+N+/OlStX8PLLL2P9+vXuDsVtysrKcOHCBeTk5ECn0+Gee+7BN77xDahUveerJSIiAgkJCdi9ezcAYNeuXUhNTXVvUC4SHR0NvV7f5jZXflb2nt8yJ7NarViwYAFefPFFxMTEuDsct4iPj8fTTz8NADAajdiwYQPmzp3r5qiU4chp7L4oODgYS5cuhUajQVNTE/70pz/57L95d55++mk8+eSTuP32290dituUlpZCpVLhs88+w9GjR3H77bdj/vz57g7LpdRqNd577z08+eSTCAwMxKuvvoq33nrL3WG5jSs/K5nAOOgvf/kLJEnCT37yE3eH4nZ1dXWYPn06Ro4ciQULFrg7HMXYHDiN3Vc1NzfjBz/4AVQqFVauXOnucFxu165dOHz4MJ577jl3h+JWdXV1sFgsyMjIgF6vx1NPPYW9e/eisbHR3aG5TENDA374wx9ixYoVyM/Px5w5c3rNGpjOuOqzkgmMg9577z0cO3YMkZGRGDBgAABgxIgROHDggJsjc636+npMnToVWq0W7777rs8OHd96Qrqc09h9ldVqxezZs3Hu3Dl88MEHvXIK6V//+hdKS0sRHx+P8PBwGI1GTJs2De+88467Q3Op1hODIyMjAQBRUVGw2Ww+vRvxVh9++CEaGxvxq1/9CsOGDcPq1auRl5eHY8eOuTs0t3DlZ6Vvftu4wP/93//h9OnTOHLkCPbv3w8A2LlzJ8aMGePmyFxr4cKFCA4OxrZt2xAUFOTucBTjyGnsvmrlypU4d+4c9uzZY//i6m1Wr16NM2fO4MiRIzhy5AhCQ0OxYcMGfPvb33Z3aC6VnJwMf39/nD17FgBQXl4OPz+/NlMIvs7Pzw8Wi8V+3WazwWq12teA9Dau/KzkWUhOYDAYEBERgYsXL7Zb0OTLjh49invuuQeHDh1qsw4oPDzcfUEpxGq1YvDgwZg9ezYWLFiAhQsXQqfT4Z///Ke7Q3OpsrIyDBkyBLt27Wqz9iMsLMxnR9/kCA8Px9atWzFx4kR3h+JyDz/8MGpqavD666/jmWeegclkwvbt290dlstUVlZi0KBBWLFiBb73ve9h7dq1eP/991FYWAh/f393h+cSkiShoKAAqampLv2s7L2fONRjW7duhdFoREpKCiIiIuwXX6RSqbB582ZkZWUhOTkZJpMJa9ascXdYLrd7927U1NRg/Pjxbf7Ne1PhLmrrtddeg81mwx133IHy8nK8+uqr7g7JpWJiYrB582b8/e9/x5AhQ3DgwAH85z//6TXJy61c+VnJERgiIiLyOhyBISIiIq/DBIaIiIi8DhMYIiIi8jpMYIiIiMjrMIEhIiIir8MEhoiIiLwOExgiIiLyOkxgiIiIyOswgSEiIiKvwwSGiIiIvA4TGCIiIvI6/x+Lg7m9ad0ALgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "\n",
    "def norm_dist_prob(theta):\n",
    "    y = norm.pdf(theta, loc=3, scale=2)\n",
    "    return y\n",
    "\n",
    "\n",
    "T = 5000\n",
    "pi = [0 for i in range(T)]\n",
    "sigma = 1\n",
    "t = 0\n",
    "while t < T - 1:\n",
    "    t = t + 1\n",
    "    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1,\n",
    "                       random_state=None)  #状态转移进行随机抽样\n",
    "    alpha = min(\n",
    "        1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))  #alpha值\n",
    "\n",
    "    u = random.uniform(0, 1)\n",
    "    if u < alpha:\n",
    "        pi[t] = pi_star[0]\n",
    "    else:\n",
    "        pi[t] = pi[t - 1]\n",
    "\n",
    "plt.scatter(pi, norm.pdf(pi, loc=3, scale=2), label='Target Distribution')\n",
    "num_bins = 50\n",
    "plt.hist(pi,\n",
    "         num_bins,\n",
    "         density=1,\n",
    "         facecolor='red',\n",
    "         alpha=0.7,\n",
    "         label='Samples Distribution')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b195b03a",
   "metadata": {},
   "source": [
    "# 二维Gibbs采样实例python实现  \n",
    "\n",
    "假设我们要采样的是一个二维正态分布 $N(\\mu, \\Sigma)$ ，其中： $\\mu=(\\mu_{1}, \\mu_{2})= (5, -1)$ , $\\Sigma = \\begin{pmatrix}\n",
    "\\sigma^{2}_{1} &   \\rho \\sigma_{1}\\sigma_{2}b\\rho \\sigma_{2}& \n",
    "\\sigma^{2}_{2}\\end{pmatrix} = \\begin{pmatrix}\n",
    " 1& 1b1 & \n",
    "4\\end{pmatrix}$;\n",
    "\n",
    "而采样过程中的需要的状态转移条件分布为：\n",
    "\n",
    "$P(x_{1}|x_{2}) = N(\\mu_{1}+ \\rho \\sigma_{1}/\\sigma_{2}(x_{2} - \\mu_{2}), (1 - \\rho^{2})\\sigma^{2}_{1})$\n",
    "\n",
    "$P(x_{2}|x_{1}) = N(\\mu_{2}+ \\rho \\sigma_{2}/\\sigma_{1}(x_{1} - \\mu_{1}), (1 - \\rho^{2})\\sigma^{2}_{2})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8845252f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGvCAYAAABB3D9ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0OklEQVR4nO3df3RU9Z3/8dcME5JJIEOE2AQSDYioqGmCbFCaUAoaVFAsbBEVXPyxPfxo14Uqu6GVjFRMLWsWu6Xo1lqrqF1oDdqA5eyeCtUEzlIScRcixpWYhB9JTsJkMJiQTO73D75Mze9JyK9P5vk45x65c+/7zufDzZCXn/uZe22WZVkCAAAwkH2gGwAAANBTBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQCSpL1798pms+mPf/xji9dtNps2b97cZj8AGAwIMgC6ZcqUKdq/f3+P60tKSuR2u3uvQQCCGkEGQLdERkbq5ptv7nF9SUmJnnrqqV5sEYBgRpABAADGIsgAaKG5uVlNTU3+pbXO5sh89NFHmjVrliIjI3X55Zdr+fLlamhokCS53W7ZbDZ961vfknRh7o3NZtOyZctaHGPbtm2aNGmSQkNDlZKSovfff7/F9tLSUs2ZM0cjR47UtGnTlJ2drUmTJunpp59u0b4zZ85o8+bNuvbaa7V69eoWx3j99dc1efJkhYeH69prr9Wbb77p32az2ZSVlaUrr7xScXFx2rdvn2644QZddtll+sMf/tC9v0wAfY4gA6CFuXPnKiQkxL90x/z583X+/Hn9/ve/13PPPaff/e53ys7OliQ9+uij2r9/v7Zs2SJJ2r9/v/bv368nn3zSX//qq6/qwQcf1D333KN33nlHCQkJuvXWW3Xo0CH/Po888ohsNpvefvttXXPNNfrxj3+sV155RYsWLWrRlnXr1unFF1/Uww8/rG9/+9v+1/Py8rR06VLdcccd2r17t+699149+OCD+uyzz1q0Y+vWrTp//rzmzZundevWKTk5WS+88EK3/j4A9D3HQDcAwODys5/9TLfccot//W/+5m8CqmtqalJ5eblWrFih2267TZJ0zTXXaPjw4ZKkuLg4xcXFqb6+XpLanWezfv16PfDAA/rpT38qSbrtttv09a9/XU8//bRycnIkXQhAO3bs0KxZsxQbG6vXXntN48ePV2xsbItj5efn67//+781cuTIFq87nU5/wBk2bJhuuOEGZWVl6S9/+YsmTJggSfrRj36kO++8U9ddd50mTZqk+++/X8eOHdO+ffsC+rsA0H8IMgBauPrqqzV16tRu1zkcDq1atUpPPvmk8vLylJKSonnz5unrX/96QPVVVVX6/PPPW3yjyW6361vf+pbeeust/2vXXnutcnNz9Y1vfENvv/22oqKi9LWvfa3N8TZt2tQmxEgXvnV1/vx5/fM//7Py8vJUUFCgpqYmnTt3zr/PuHHjJF24zPTVPwMYfLi0BKDXbN68Wfv379eMGTOUl5en5ORk/eIXvwio1rKsdl+32WwttiUnJ+uVV16Ry+VSVlaWXnnlFdntbf8pS0lJafd4L774ombMmCGPx6Ply5fr6NGjuuKKKwJqI4DBhyADoFeUl5fre9/7niZPnqwf/OAH2r17txYtWqQXX3yxxX5hYWGSpC+//LLF65dffrni4+P1pz/9yf9ac3Oz3nvvPf8IUWFhod555x1VV1fr6NGjOnXqlO6+++5utfOFF17Q4sWL9ctf/lIPPvigRowYoZqamp50GcAgwKUlAL1i1KhReuONN3T+/HktXrxYtbW1ys/P10033dRiv+uuu04jRozQs88+qxkzZuh//ud/dO+99yomJkZPPfWUHn30UY0bN07f+ta39PLLL+vjjz/WSy+9JOnC/Jbq6mr94he/0E033aRz584pNjZWY8eODbido0eP1v79+7Vnzx6dPn1aGzdu1NmzZ9v9hhaAwY8RGQC9YsSIEdq1a5c++eQTffvb39ZDDz2km266yf8tpYtcLpfefPNNvfHGG5ozZ06LS08PPfSQXn75Zf3ud7/TvHnz9H//93/6z//8T/9lookTJ2rq1KnKysrS7bffrqlTp2rcuHGaMWNGmxGejvzbv/2bYmJidM8998jtdmvlypWaOnVqm695AzCDzerowjQADDI/+tGPtGfPHj399NMaOXKkGhsb9ec//1nr16/X4cOHlZiYONBNBNDPuLQEwBhLly7V4cOH9eCDD+rMmTMKDQ3Vtddeq5/97GeEGCBIMSIDAACMxRwZAABgLIIMAAAwFkEGAAAYiyADAACMZfS3lpqbm3Xy5EmNHDmS56AAAGAIy7J09uxZjR07tt1HjHSH0UHm5MmTio+PH+hmAACAHigrK1NcXNwlHcPoIHPxybZlZWWKjIwc4NYAAIBAeL1excfHt/uE+u4yOshcvJwUGRlJkAEAwDC9MS2Eyb4AAMBYBBkAAGAsggwAADCW0XNkAAAwmc/nU2Nj40A3o9cNGzZMDoejX26NQpABAGAAfPHFFyovL9dQfXZzeHi4YmNjNXz48D59H4IMAAD9zOfzqby8XOHh4YqOjh5SN3W1LEvnz59XVVWVjh8/rquvvvqSb3rXGYIMAAD9rLGxUZZlKTo6Wk6nc6Cb0+ucTqdCQkL0+eef6/z58woLC+uz92KyLwAAA2QojcS01pejMC3ep1/eBQAAoA9waQkAgEHCvdfdv+83s3/fry8wIgMAAIxFkAEAAMYiyAAAAGMRZAAAQECqq6s1ZswY/fGPf5Qk/eu//quSk5Pl8/kGrE1M9gUA9JvWk1mHwmTTYDJ69Gg988wzysjI0M0336ysrCzt3LlTw4YNG7A2MSIDAAAC9uijjyokJETp6emaO3eupk+fPqDt6VaQKSgoUFJSkpxOp9LT01VZWRlw7auvviqbzaaSkhL/a3v27NGkSZMUERGhxYsX69y5c91pDgAA6Gd2u13f/e53dfDgQS1fvnygmxN4kGlubtbChQs1b948FRcXy+l0avXq1QHVejwePfHEE21eW7RokdasWaMjR46opKREGzdu7F7rAQBAv6qvr9dPfvITzZkzR+vXrx/o5gQeZPbt26eamhq53W7FxcUpMzNTOTk5qqur67J23bp1bYaecnJyFBcXp+XLlyshIUFr167Vtm3but8DAADQb3784x/rqquu0jvvvKPPPvtMb7zxxoC2J+DJvnl5eUpJSZHDcaEkKSlJPp9PBQUFSktL67Du0KFD2rFjh/Lz87Vz584Wx/tquJk2bZpKS0tVVlam+Pj4do/V0NCghoYG/7rX6w20+QAADHqDffLz0aNH9fzzz+svf/mLhg8frs2bN+vhhx/WHXfcoaioqAFpU8BB5vTp0xozZox/3W63KyoqShUVFR3WNDc3a8WKFXr66acVHR3d5ng33nijf3306NGSpIqKig6DTFZWlp566qlAmwwAGOT4FpNZJk+erC+++MK/Pnfu3E5zQH/o1mRfy7LarHf25M4XX3xRNptNf//3f9/l8S7+ubPjZWRkqLa21r+UlZV1p/kAAGCICXhEJjY2VkVFRf51n88nj8ejmJiYDmu2b9+ujz76SJdddpk/qCQmJmr37t2KjY1VVVWVf9/q6mpJ6vR4oaGhCg0NDbTJAABgiAt4RCYtLU0HDx5UU1OTJKmwsFAOh0PJyckd1rz55psqKirShx9+qPfff1+StHv3bk2dOlVpaWnKz8/373vgwAElJCRo3LhxPe0LAAAIMgEHmdTUVEVHRyszM1Pl5eXasGGDFixYoPDwcHm9XjU2NrapiYmJUUJCghISEnTFFVdIkuLi4hQWFqb58+fr1KlT2rp1q0pKSrRp0yYtWbKk93oGAACGvICDjN1u144dO5Sbm6uJEyeqvr5e2dnZki5cLtq1a1e33tjlcmn79u3Kzs7W5MmTNX78eGVkZHSv9QAAIKh161lLU6ZM0eHDh9u8/tW79XZk1KhRbSYLp6enq7i4uDtNAAAYpPW3koDexrOWAACAsXj6NQAAg4XbPbTfrw8wIgMAAIxFkAEAAMYiyAAAAGMRZAAAQEBee+013XDDDf71L774QmFhYfr4448HrE0EGQAAEJD58+fr008/1bFjxyRJ7777rq655hpde+21A9YmggwAAAhIZGSkbr/9dv3+97+XJO3cuVP33nvvgLaJIAMAAAK2aNEivfXWW2psbNTu3bsJMgAAwBx33323jh49ql//+te66qqrdNVVVw1oe7ghHgBg0Gj9SAP3THe7+2HgjBgxQnfccYfWrl2rH/7whwPdHIIMAACDhiF32r333nv11ltvadGiRQPdFC4tAQCAwJWUlOjLL7/U9OnTdeWVVw50cxiRAQAAgbv77rtVWVmpHTt2DHRTJBFkAABAN3z00UcD3YQWuLQEAACMxYgMAKDXtP7WETpnWdZAN6HP9FffGJEBAKCfDRs2TJJ0/vz5AW5J3zl37pwkKSQkpE/fhxEZAAD6mcPhUHh4uKqqqhQSEiK7feiMK1iWpXPnzqmyslKjRo3yh7a+QpABAKCf2Ww2xcbG6vjx4/r8888Hujl9YtSoUYqJienz9yHIAAAwAIYPH66rr756SF5eCgkJ6fORmIsIMgAADBC73a6wsLCBbobRhs5FOQAAEHQIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsQgyAADAWN0KMgUFBUpKSpLT6VR6eroqKyu7rMnLy9ONN96okSNHau7cuS1qYmJiZLPZ/MvUqVO73wMAwJDl3utuswBfFXCQaW5u1sKFCzVv3jwVFxfL6XRq9erVndbU19frO9/5jp544gkVFRVJkjIyMvzbPR6PioqKdObMGZ05c0Z79+7tWS8AAEBQCvihkfv27VNNTY3cbrccDocyMzOVmpqquro6RUREtFtz/PhxTZ06VQ8++KAkae7cuXrttdckSQ0NDWpoaNCECRM0fPjwXugKAAAINgGPyOTl5SklJUUOx4Xsk5SUJJ/Pp4KCgg5rrrvuOr3zzjuSpLNnz2rnzp1KTk6WdGE0JiwsTPfdd5+cTqdmzJihEydOdNqGhoYGeb3eFgsAAAheAQeZ06dPa8yYMX8ttNsVFRWlioqKLmtLS0s1atQolZSUaOPGjZIuBJn6+nrNnj1bR48eld1u1+OPP97pcbKysuRyufxLfHx8oM0HAABDULcm+1qW1WbdZrN1WTd27Fh98MEHcrlcevLJJyVJEyZM0MmTJ7Vy5UqNHz9eK1as0HvvvdfpcTIyMlRbW+tfysrKutN8AAAwxAQcZGJjY1VVVeVf9/l88ng8iomJ6bCmrKxMRUVFcjgcuuWWW7Ru3Tq9/vrrkqSQkBDFxsb6942KiuryUlFoaKgiIyNbLAAAIHgFHGTS0tJ08OBBNTU1SZIKCwvlcDj8c17a8+677+qRRx7xr9tsNg0bNkyStGXLFs2aNcu/rbS0VAkJCd1tPwAACGIBB5nU1FRFR0crMzNT5eXl2rBhgxYsWKDw8HB5vV41Nja2qZk1a5YKCwuVk5Ojzz//XFu2bNGdd94pSZo5c6by8/OVm5urY8eO6fnnn9eyZct6rWMAAGDoC/jr13a7XTt27NBDDz2k5557TjNmzNCvfvUrSVJiYqI2b96se+65p0XNxIkT9eqrr2rt2rWqrKzU7bffrs2bN0uSrr/+er3wwgtatWqVvvzySz3wwANd3pcGADB4cHM6DAY2q/UMXoN4vV65XC7V1tYyXwYA+tlABRn3zIF5X/Se3vz9zbOWAACAsQgyAADAWAQZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEcA90AAAC6w73X3XJ9prvd/RAcGJEBAADGIsgAAABjEWQAAICxCDIAAMBYTPYFAASk9SRbYDBgRAYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGKtbQaagoEBJSUlyOp1KT09XZWVllzV5eXm68cYbNXLkSM2dO7dFzZ49ezRp0iRFRERo8eLFOnfuXPd7AAAAglbAQaa5uVkLFy7UvHnzVFxcLKfTqdWrV3daU19fr+985zt64oknVFRUJEnKyMiQJHk8Hi1atEhr1qzRkSNHVFJSoo0bN15CVwAAQLAJOMjs27dPNTU1crvdiouLU2ZmpnJyclRXV9dhzfHjxzV16lQ9+OCDiouL09y5c3X06FFJUk5OjuLi4rR8+XIlJCRo7dq12rZt26X3CAAABI2Ag0xeXp5SUlLkcFx4zmRSUpJ8Pp8KCgo6rLnuuuv0zjvvSJLOnj2rnTt3Kjk52X+86dOn+/edNm2aSktLVVZW1qOOAACA4BNwkDl9+rTGjBnz10K7XVFRUaqoqOiytrS0VKNGjWpx+aj18UaPHi1JnR6voaFBXq+3xQIAAIJXtyb7WpbVZt1ms3VZN3bsWH3wwQdyuVx68skn2z3exT93drysrCy5XC7/Eh8f353mAwCAISbgIBMbG6uqqir/us/nk8fjUUxMTIc1ZWVlKioqksPh0C233KJ169bp9ddfb/d41dXVktTp8TIyMlRbW+tfuAwFAEBwCzjIpKWl6eDBg2pqapIkFRYWyuFw+Oe8tOfdd9/VI4884l+32WwaNmyY/3j5+fn+bQcOHFBCQoLGjRvX4fFCQ0MVGRnZYgEAAMEr4CCTmpqq6OhoZWZmqry8XBs2bNCCBQsUHh4ur9erxsbGNjWzZs1SYWGhcnJy9Pnnn2vLli268847JUnz58/XqVOntHXrVpWUlGjTpk1asmRJ7/UMAAAMeQEHGbvdrh07dig3N1cTJ05UfX29srOzJUmJiYnatWtXm5qJEyfq1Vdf1dq1a5WYmKjLLrtMmzdvliS5XC5t375d2dnZmjx5ssaPH++/xwwAAEAgbFbrGbwG8Xq9crlcqq2t5TITAPQx9173QDehXe6Z7oFuArqpN39/86wlAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsx0A3AACAS9H6jsPc6Te4MCIDAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMJZjoBsAABic3HvdA90EoEuMyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjNWtIFNQUKCkpCQ5nU6lp6ersrKyy5pDhw5pypQpioyM1Pz581VTU+PfFhMTI5vN5l+mTp3a/R4AAICgFXCQaW5u1sKFCzVv3jwVFxfL6XRq9erVXdbcf//9Sk9P10cffaTKykqtX7/ev93j8aioqEhnzpzRmTNntHfv3h53BAAABJ+A7+y7b98+1dTUyO12y+FwKDMzU6mpqaqrq1NERES7NZ999pk++eQTPfnkk4qIiNDSpUv14osvSpIaGhrU0NCgCRMmaPjw4b3TGwAAEFQCHpHJy8tTSkqKHI4L2ScpKUk+n08FBQUd1oSHh2vz5s3+oFNdXS2n0ynpwmhMWFiY7rvvPjmdTs2YMUMnTpzotA0NDQ3yer0tFgAAELwCDjKnT5/WmDFj/lpotysqKkoVFRUd1owdO1aPPfaYJKm2tlYvvfSSli5dKulCkKmvr9fs2bN19OhR2e12Pf744522ISsrSy6Xy7/Ex8cH2nwAADAEdWuyr2VZbdZtNluXdXV1dbrrrruUnJysFStWSJImTJigkydPauXKlRo/frxWrFih9957r9PjZGRkqLa21r+UlZV1p/kAAGCICXiOTGxsrIqKivzrPp9PHo9HMTExndadO3dOc+bM0YgRI/Tb3/5WdvuF7BQSEqLY2Fj/flFRUV1eKgoNDVVoaGigTQYABKHWT+12z3S3ux+GhoBHZNLS0nTw4EE1NTVJkgoLC+VwOJScnNxp3cqVKxUREaG3335bYWFh/te3bNmiWbNm+ddLS0uVkJDQzeYDAIBgFnCQSU1NVXR0tDIzM1VeXq4NGzZowYIFCg8Pl9frVWNjY5uaw4cPa+fOnfr5z3+uL7/8Uh6PRx6PR5I0c+ZM5efnKzc3V8eOHdPzzz+vZcuW9Va/AABAEAg4yNjtdu3YsUO5ubmaOHGi6uvrlZ2dLUlKTEzUrl272tTs3LlTtbW1mjRpkqKiovyLJF1//fV64YUXtGrVKqWlpenWW2/t8r40AAAAX2WzWs/gNYjX65XL5VJtba0iIyMHujkAMKS0nmtiKubIDD69+fubZy0BAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYK+KGRAIChbajcAA/BhREZAABgLIIMAAAwFkEGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsboVZAoKCpSUlCSn06n09HRVVlZ2WXPo0CFNmTJFkZGRmj9/vmpqavzb9uzZo0mTJikiIkKLFy/WuXPnut8DAAAQtAIOMs3NzVq4cKHmzZun4uJiOZ1OrV69usua+++/X+np6froo49UWVmp9evXS5I8Ho8WLVqkNWvW6MiRIyopKdHGjRsvrTcAACCoBBxk9u3bp5qaGrndbsXFxSkzM1M5OTmqq6vrsOazzz7TJ598oieffFIJCQlaunSp3n//fUlSTk6O4uLitHz5ciUkJGjt2rXatm3bpfcIAAAEjYCDTF5enlJSUuRwOCRJSUlJ8vl8Kigo6LAmPDxcmzdvVkREhCSpurpaTqfTf7zp06f79502bZpKS0tVVlbW4fEaGhrk9XpbLAAAIHg5At3x9OnTGjNmjH/dbrcrKipKFRUVHdaMHTtWjz32mCSptrZWL730ktauXes/3o033ujfd/To0ZKkiooKxcfHt3u8rKwsPfXUU4E2GcClcrsvbTswCLj3uluuz3S3ux/MFHCQkSTLstqs22y2Luvq6up01113KTk5WStWrGj3eBf/3NnxMjIytGbNGv+61+vtMPQACABBBIDhAg4ysbGxKioq8q/7fD55PB7FxMR0Wnfu3DnNmTNHI0aM0G9/+1vZ7Xb/8aqqqvz7VVdXS1KnxwsNDVVoaGigTQYAAENcwHNk0tLSdPDgQTU1NUmSCgsL5XA4lJyc3GndypUrFRERobffflthYWEtjpefn+9fP3DggBISEjRu3Lju9gEAAASpgINMamqqoqOjlZmZqfLycm3YsEELFixQeHi4vF6vGhsb29QcPnxYO3fu1M9//nN9+eWX8ng88ng8kqT58+fr1KlT2rp1q0pKSrRp0yYtWbKk1zoGAACGvoAvLdntdu3YsUMPPfSQnnvuOc2YMUO/+tWvJEmJiYnavHmz7rnnnhY1O3fuVG1trSZNmtTidcuy5HK5tH37dq1atUo/+MEPdPfddysjI+PSewSg/zAZGMAAs1mtZ/AaxOv1yuVyqba2VpGRkQPdHMA8Ax00Bvr90ULrb/cMVXxraeD15u9vnrUEAACMRZABAADG6tZ9ZAAAQ0ewXErC0EaQAYYy5qAAGOK4tAQAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxuIRBQAGTlePUOARCwC6wIgMAAAwFkEGAAAYiyADAACMRZABAADGYrIvYDImwwIIcozIAAAAYxFkAACAsbi0BAAIKu697pbrM93t7gczMCIDAACMxYgMMJgxmRcAOsWIDAAAMBZBBgAAGItLSwAGr0AurXH5DQhqjMgAAABjEWQAAICxCDIAAMBYBBkAAGCsbgWZgoICJSUlyel0Kj09XZWVlQHVffrpp5o8ebI8Hk+L12NiYmSz2fzL1KlTu9McAAAQ5AIOMs3NzVq4cKHmzZun4uJiOZ1OrV69usu6OXPm6Oqrr1ZRUVGbbR6PR0VFRTpz5ozOnDmjvXv3dqvxAAAguAX89et9+/appqZGbrdbDodDmZmZSk1NVV1dnSIiIjqs+81vfqNPP/1UaWlpLV5vaGhQQ0ODJkyYoOHDh/e8BwCAgLR+xhAwFAQ8IpOXl6eUlBQ5HBeyT1JSknw+nwoKCjqti4mJUVxcXJvXPR6PwsLCdN9998npdGrGjBk6ceJEp8dqaGiQ1+ttsQAAgOAVcJA5ffq0xowZ89dCu11RUVGqqKjo0Rt7PB7V19dr9uzZOnr0qOx2ux5//PFOa7KysuRyufxLfHx8j94bAAAMDd2a7GtZVpt1m83WozeeMGGCTp48qZUrV2r8+PFasWKF3nvvvU5rMjIyVFtb61/Kysp69N4AAGBoCHiOTGxsbIsJuz6fTx6PRzExMT1645CQEMXGxvrXo6KiurxUFBoaqtDQ0B69HwAAGHoCHpFJS0vTwYMH1dTUJEkqLCyUw+FQcnJyj954y5YtmjVrln+9tLRUCQkJPToWAAAITgEHmdTUVEVHRyszM1Pl5eXasGGDFixYoPDwcHm9XjU2NnbrjWfOnKn8/Hzl5ubq2LFjev7557Vs2bLuth8AAASxgIOM3W7Xjh07lJubq4kTJ6q+vl7Z2dmSpMTERO3atatbb3z99dfrhRde0KpVq5SWlqZbb701oPvSAAAAXGSzWs/gNYjX65XL5VJtba0iIyMHujlA73O7B7oF5uPv0I/7yLTPPdM90E0IOr35+5tnLQEAAGMRZAAAgLEIMgAAwFgB30cGQB9g/gYAXBJGZAAAgLEIMgAAwFhcWgIABLXWX0vn69hmYUQGAAAYiyADAACMRZABAADGIsgAAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIzFIwqAvsTTrQGgTxFkAAxtXYVJwiZgNC4tAQAAYxFkAACAsQgyAADAWMyRAYAhyL3XPdBNAPoFIzIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgAAwFjdCjIFBQVKSkqS0+lUenq6KisrA6r79NNPNXnyZHk8nhav79mzR5MmTVJERIQWL16sc+fOdac5AAAgyAUcZJqbm7Vw4ULNmzdPxcXFcjqdWr16dZd1c+bM0dVXX62ioqIWr3s8Hi1atEhr1qzRkSNHVFJSoo0bN3a/BwAAIGgFHGT27dunmpoaud1uxcXFKTMzUzk5Oaqrq+u07je/+Y3ef//9Nq/n5OQoLi5Oy5cvV0JCgtauXatt27Z1vwcAACBoBfzQyLy8PKWkpMjhuFCSlJQkn8+ngoICpaWldVgXExOj+vr6do83ffp0//q0adNUWlqqsrIyxcfHd6cPwMBxuwe6BQAQ1AIekTl9+rTGjBnz10K7XVFRUaqoqOjRG7c+3ujRoyWp0+M1NDTI6/W2WAAAQPAKeERGkizLarNus9l6/OZfPd7FP3d2vKysLD311FM9fj8AaKOrUTVG3YBBLeAgExsb22LCrs/nk8fjUUxMTI/eODY2VlVVVf716upqSer0eBkZGVqzZo1/3ev1chkKANCr3HvdbV+b2fY1DA4BX1pKS0vTwYMH1dTUJEkqLCyUw+FQcnJyj944LS1N+fn5/vUDBw4oISFB48aN67AmNDRUkZGRLRYAABC8Ag4yqampio6OVmZmpsrLy7VhwwYtWLBA4eHh8nq9amxs7NYbz58/X6dOndLWrVtVUlKiTZs2acmSJd3uAAAACF4BBxm73a4dO3YoNzdXEydOVH19vbKzsyVJiYmJ2rVrV7fe2OVyafv27crOztbkyZM1fvx4ZWRkdK/1AAAgqHVrsu+UKVN0+PDhNq+XlJR0WpeQkNBmorAkpaenq7i4uDtNAAAA8ONZSwAAwFjdGpEBAAxO7X3TBggGjMgAAABjEWQAAICxCDIAAMBYBBkAAGAsJvsCQGd4FhMwqDEiAwAAjMWIDNAR/k8bAAY9RmQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGNxQzwAALrg3utuuT7T3e5+6H+MyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCwm+wLApejqKek8RR3oUwQZBC9+wcBgrb9FAwQrLi0BAABjEWQAAICxCDIAAMBYBBkAAGAsggwAADAWQQYAABiLIAMAAIxFkAEAAMYiyAAAAGN1686+BQUFevjhh3Xs2DGlpaVp27Ztuvzyyzut2bNnj77//e/rxIkTuuuuu/Tyyy8rPDxckhQTE6OKigr/vjfddJP+8pe/9KAbADBI8QgDoE8FHGSam5u1cOFCPfDAA8rNzdWqVau0evVqvf766x3WeDweLVq0SM8++6xuv/12LV68WBs3btTGjRv924uKihQTE3OhMQ6emIBexC8IABjyAk4O+/btU01NjdxutxwOhzIzM5Wamqq6ujpFRES0W5OTk6O4uDgtX75ckrR27VqtXr1aGzduVENDgxoaGjRhwgQNHz68d3oDAEA/aP2sK/dMd7v7oe8FPEcmLy9PKSkp/lGTpKQk+Xw+FRQUdFozffp0//q0adNUWlqqsrIyeTwehYWF6b777pPT6dSMGTN04sSJTtvQ0NAgr9fbYgEAAMEr4CBz+vRpjRkz5q+FdruioqJazHHpqmb06NGSpIqKCnk8HtXX12v27Nk6evSo7Ha7Hn/88U7bkJWVJZfL5V/i4+MDbT4AABiCuvWtJcuy2qzbbLaAay7+2WazacKECTp58qRWrlyp8ePHa8WKFXrvvfc6PVZGRoZqa2v9S1lZWXeaDwAAhpiA58jExsaqqKjIv+7z+eTxePwTdTuqqaqq8q9XV1dLuvBtpZCQEMXGxvq3RUVFdXmpKDQ0VKGhoYE2GQAADHEBj8ikpaXp4MGDampqkiQVFhbK4XAoOTm505r8/Hz/+oEDB5SQkKBx48Zpy5YtmjVrln9baWmpEhISetAFAAAQrAIekUlNTVV0dLQyMzO1YsUKbdiwQQsWLFB4eLi8Xq+cTqdCQkJa1MyfP1//8A//oK1bt+qOO+7Qpk2btGTJEknSzJkz9YMf/EC5ubm6+uqr9fzzz2vZsmW92jkAGCpaf0sGwAUBj8jY7Xbt2LFDubm5mjhxourr65WdnS1JSkxM1K5du9rUuFwubd++XdnZ2Zo8ebLGjx+vjIwMSdL111+vF154QatWrVJaWppuvfVWrV69upe6BQAAgkG37kA3ZcoUHT58uM3rJSUlHdakp6eruLi43W3Lli1jFAYAAPQYz1oCAADG4pkAADCQeBYTcEkIMjAX/8ADQNDj0hIAADAWQQYAABiLS0sAAFwinoY9cBiRAQAAxmJEBgAGob0ley/8lzv6Ap1iRAYAABiLIAMAAIxFkAEAAMYiyAAAAGMx2ReDF3fuBTTzlb2dbt+7bGa/tAMYrBiRAQAAxiLIAAAAYxFkAACAsZgjAwCDwMUb4AHoHoIMAAC9jGcv9R8uLQEAAGMxIoOBw9ergUvG17MR7BiRAQAAxmJEBgCAPsacmb7DiAwAADAWQQYAABiLS0voO0zmBQYck4Ex1BFkAKCfcfM7oPdwaQkAABiLIAMAAIxlsyzLGuhG9JTX65XL5VJtba0iIyMHujnBhzkwQEBMvpTEHJr+E0xfye7N39+MyAAAAGMx2RcdY8QFADDIEWQAoJeZfCmpta6+vi1x+am3cPffnulWkCkoKNDDDz+sY8eOKS0tTdu2bdPll1/eac2ePXv0/e9/XydOnNBdd92ll19+WeHh4V1uQz9gxAUABi2CTWACnuzb3Nysq666Sg888ICWL1+uVatWacSIEXr99dc7rPF4PLryyiv17LPP6vbbb9fixYs1e/Zsbdy4sdNtgWKy7yUiyAC9YiiNwPQFRmz6hsnBpjd/fwccZN577z3dc889qq6ulsPhUEFBgVJTU1VVVaWIiIh2a37961/rX/7lX3TkyBFJ0ltvvaXVq1fr888/73RboAgyXSCoAH2C4NK7CDo9Q5C5IOBLS3l5eUpJSZHDcaEkKSlJPp9PBQUFSktL67Bm+vTp/vVp06aptLRUZWVlnW6Lj49v93gNDQ1qaGjwr9fW1kq68BcSlLKyBroFQFCqO9800E0YUv7m3/+r0+0fPND+75hgl7Ero/PtaZ1vH0gXf2/3xh1gAg4yp0+f1pgxY/zrdrtdUVFRqqio6LTmxhtv9K+PHj1aklRRUdHpto6CTFZWlp566qk2r3e0PwBgCHgzb6BbYKSf6CcD3YQunT17Vi6X65KO0a3Jvq2Tk2VZstlsAddc/PPFms62tScjI0Nr1qzxrzc3N6umpkajR4/ush1DhdfrVXx8vMrKyoLmchp9Do4+S8HZ72DssxSc/Q7GPkvt99uyLJ09e1Zjx4695OMHHGRiY2NVVFTkX/f5fPJ4PIqJiem0pqqqyr9eXV0tSYqJiel0W0dCQ0MVGhra4rVRo0YF2oUhJTIyMqg+CBJ9DibB2O9g7LMUnP0Oxj5Lbft9qSMxFwV8Z9+0tDQdPHhQTU0Xrg0XFhbK4XAoOTm505r8/Hz/+oEDB5SQkKBx48Z1ug0AACAQAQeZ1NRURUdHKzMzU+Xl5dqwYYMWLFig8PBweb1eNTY2tqmZP3++Tp06pa1bt6qkpESbNm3SkiVLutwGAAAQiICDjN1u144dO5Sbm6uJEyeqvr5e2dnZkqTExETt2rWrTY3L5dL27duVnZ2tyZMna/z48crIyOhyGzoWGhqqzMzMNpfYhjL6HDyCsd/B2GcpOPsdjH2W+r7fRj/9GgAABDeefg0AAIxFkAEAAMYiyAAAAGMRZAAAgLEIMgbYu3evbDZbiyUhIaHLuoaGhjZ1f/u3f9v3De5FMTExLdo/derUgOpee+01xcfHa+TIkfr+978vn8/Xxy3tPYcOHdKUKVMUGRmp+fPnq6ampssaU851QUGBkpKS5HQ6lZ6ersrKyi5r9uzZo0mTJikiIkKLFy/WuXPn+qGlveuzzz7TN7/5TY0cOVIzZ84M6OG4Pf3ZHyxuvvnmFu3/6iNuOmL6uS4pKWnzOQzkrvMmnutTp07pm9/8pj788EP/az35fEu9cN4tDHqNjY3WmTNn/MvTTz9t3XzzzV3WnTp1ygoJCWlRW1dX1w8t7j2hoaFWUVGRv/1nz57tsubYsWNWWFiYtXPnTquoqMiaMGGC9e///u/90NpL5/P5rEmTJln/9E//ZB0/fty6+eabrVWrVnVZZ8K59vl8VkJCgvXDH/7QKisrs+6++27r/vvv77TmzJkzVmRkpLV161br+PHj1rRp06x169b1U4t7z+zZs60lS5ZYJSUl1oIFC6y5c+d2WdOTn/3B5JprrrH27Nnjb39tbW2n+w+Fc+3z+Vp8Bl999VUrJiamyzrTzvV3v/tdS5IlySosLLQsq2efb8vqnfNOkDHQrbfearnd7i73+/jjj62xY8f2Q4v6Rn19vSXJamho6Fbd+vXrrTvvvNO//txzz1lpaWm93bw+UVxcbEmyvvjiC8uyLGvLli1WYmJil3UmnOs//elPVmRkpNXY2GhZlmUdOnTIcjqd/r625+WXX7YmT57sX//9739vXXHFFX3e1t7U0NBg2Ww268iRI5ZlWdauXbusyMjITmt6+rM/mMTExFiffPJJwPsPhXPd2qOPPmotW7as031MPNdVVVXW8ePHWwSZnny+Lat3zjuXlgzzxRdf6M9//rPuvPPOLvf1eDwaNmyYZs6cqfDwcM2fP1+1tbX90Mre4fF4FBYWpvvuu09Op1MzZszQiRMnuqzLy8vT9OnT/evTpk3T/v37e+Vx8X0tPDxcmzdvVkREhKQLzyBzOp1d1plwrvPy8pSSkiKH48Ij3pKSkuTz+VRQUNBpTetzWVpaqrKysj5vb29pbGzUT3/6U40fP15SYOe0pz/7g4nH49H69evldDqVlJSkI0eOdLr/UDjXrb377rtd/ltt4rkeM2ZMm+kNPfl8X6y71PNOkDHMf/3Xf8nlcgV0DdXj8aiiokKrVq3S4cOHVVxcrJ/8ZPA/1v0ij8ej+vp6zZ49W0ePHpXdbtfjjz/eZd3p06dbXI8fPXq0mpqa/A8mHczGjh2rxx57TJJUW1url156SUuXLu2yzoRz3fq82O12RUVFqaKiIuCa0aNHS1KnNYNNRESEHn/8cTmdTjU2NupnP/tZl+e0pz/7g8X58+dVX1+viRMn6ujRo5o8ebIeffTRTmuGwrn+qg8//FAVFRW67bbbOt3P9HN9UU8+3+3V9eS8B/z0awwOu3fv1u233x7QBLJvfvObKi8vV3R0tCTp7/7u75STk9PXTew1EyZM0MmTJxUbGytJWrFihf+XfFe+Ovpy8c+B/J0NFnV1dbrrrruUnJysFStWdLm/Kee69aiYZVldnhfTz+VFTU1NeuCBB2S327Vhw4ZO972Un/3BICQkROXl5f6HAD/22GO6+eab9eWXX3Y6GjVUzrV04d/qW265RaNGjep0P9PP9Vf15PPduq4n550RGcO8++67uuOOOwLaNywszP+LTZKioqLk9Xr7qmm9LiQkxP/hlgJvf2xsrKqqqvzr1dXVCgkJ0WWXXdYn7ext586d05w5cxQeHq7f/va3stu7/piacK5bnxefzyePx6OYmJiAay6OqnVWMxg1Nzdr8eLF+vTTT/Xuu+92eWmppz/7g4XNZvOHGOlC+yXp7NmzHdYMlXN90e7duwP6t9r0c31RTz7f7dX15LwTZAxy+PBhnTx5UnPmzAlo/yeeeEIPP/ywf720tDSgr20PFlu2bNGsWbP864G2Py0tTfn5+f71AwcO6Bvf+IYx/2e3cuVKRURE6O2331ZYWFhANSac67S0NB08eFBNTU2SpMLCQjkcDiUnJ3da0/pcJiQktPglaYINGzbo008/1Z/+9KeAAnVPf/YHi9zcXF111VX+9dLSUoWHh7cI260NlXMtSWfOnNGBAwcCmsto+rm+qCef74t1l3zeuzU1GAPqmWeesW655ZY2r585c8Zqampq8/ru3butUaNGWR988IF16NAha+zYsdZ//Md/9EdTe8X//u//WqGhodYf/vAH6+OPP7ZuuOEG69lnn/Vv76jfxcXFLb5+PX78eOull17qz6b32Icffmi5XC7rk08+afE1zotMPtc+n8+aMGGCtW7dOqusrMy66667rCVLlliWZVm1tbXW+fPn29R4PB7L5XJZv/jFL6zjx49bKSkp1o9+9KP+bvolOXXqlBUZGWnl5+e3OKc+n6/Dfnf1sz/YVVZWWiNGjLB++ctfWp999pk1a9Ysa8WKFZZlDe1zfdEbb7zR5luEQ/Fcq9XXrzv6fFtW3553goxBUlNTrR//+MdtXv/qD1NrzzzzjPW1r33NGjt2rPXMM8/0cQt7369//WvriiuusKKjo61//Md/bPFB6Kzfr732mjVu3DhrxIgR1ve+9z3L5/P1U4svjdvt9t+f4avLRaaf60OHDlmJiYlWaGioddttt1mVlZWWZVnWlVdeaeXk5LRbs2fPHmvixImW0+m07r333kF3f5yuvPLKK+2e0+PHj3fa785+9k2wa9cua9KkSdaoUaOspUuX+u+NMpTP9UVLliyxHnnkkRavDcVz3frfo44+35bVt+fd9v8bAwAAYBzmyAAAAGMRZAAAgLEIMgAAwFgEGQAAYCyCDAAAMBZBBgAAGIsgAwAAjEWQAQAAxiLIAAAAYxFkAACAsf4fWSbGmnEES4IAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from scipy.stats import multivariate_normal\n",
    "\n",
    "samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])\n",
    "\n",
    "def p_ygivenx(x, m1, m2, s1, s2):\n",
    "    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))\n",
    "\n",
    "def p_xgiveny(y, m1, m2, s1, s2):\n",
    "    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))\n",
    "\n",
    "N = 5000\n",
    "K = 20\n",
    "x_res = []\n",
    "y_res = []\n",
    "z_res = []\n",
    "m1 = 5\n",
    "m2 = -1\n",
    "s1 = 1\n",
    "s2 = 2\n",
    "\n",
    "rho = 0.5\n",
    "y = m2\n",
    "\n",
    "for i in range(N):\n",
    "    for j in range(K):\n",
    "        x = p_xgiveny(y, m1, m2, s1, s2)   #y给定得到x的采样\n",
    "        y = p_ygivenx(x, m1, m2, s1, s2)   #x给定得到y的采样\n",
    "        z = samplesource.pdf([x,y])\n",
    "        x_res.append(x)\n",
    "        y_res.append(y)\n",
    "        z_res.append(z)\n",
    "\n",
    "num_bins = 50\n",
    "plt.hist(x_res, num_bins,density=1, facecolor='green', alpha=0.5,label='x')\n",
    "plt.hist(y_res, num_bins, density=1, facecolor='red', alpha=0.5,label='y')\n",
    "plt.title('Histogram')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
